The Effects of Climate Change on the Breeding Success and Timing of African White-backed Vultures

Author

Mieke Deyzel

#remotes::install_github("ropensci/rnoaa")
#remotes::install_github("ropensci/chirps")
#install.packages("GSODR")
#install.packages("AICcmodavg")

Load Libraries

library(tidyverse)
── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
✔ dplyr     1.1.4     ✔ readr     2.1.5
✔ forcats   1.0.0     ✔ stringr   1.5.1
✔ ggplot2   4.0.1     ✔ tibble    3.2.1
✔ lubridate 1.9.3     ✔ tidyr     1.3.1
✔ purrr     1.0.2     
── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
✖ dplyr::filter() masks stats::filter()
✖ dplyr::lag()    masks stats::lag()
ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
library(lubridate)
library(janitor)

Attaching package: 'janitor'

The following objects are masked from 'package:stats':

    chisq.test, fisher.test
library(MASS)

Attaching package: 'MASS'

The following object is masked from 'package:dplyr':

    select
library(patchwork)

Attaching package: 'patchwork'

The following object is masked from 'package:MASS':

    area
library(stringr)
library(tibble) 
library(dplyr)
library(knitr)
library(ggplot2)
library(GGally)
library(tidyr)
library(chirps)
library(GSODR)
library(tibble)
library(knitr)
library(AICcmodavg)
library(lme4)
Loading required package: Matrix

Attaching package: 'Matrix'

The following objects are masked from 'package:tidyr':

    expand, pack, unpack


Attaching package: 'lme4'

The following object is masked from 'package:AICcmodavg':

    checkConv

Overdispersion Check

check_overdispersion <- function(model) {
  rp <- residuals(model, type = "pearson")
  disp <- sum(rp^2) / model$df.residual
  disp
} 
vultures    <- read.csv("~/Desktop/Vulture Project/Data/93_24 cleaned.csv")
# WeatherData <- read.csv("~/Desktop/Vulture Project/Data/Weather.csv")

Weather Data

library(rnoaa)
Registered S3 method overwritten by 'hoardr':
  method           from
  print.cache_info httr
The rnoaa package will soon be retired and archived because the underlying APIs have changed dramatically. The package currently works but does not pull the most recent data in all cases. A noaaWeather package is planned as a replacement but the functions will not be interchangeable.
library(tidyverse)
library(lubridate)

 
kim_station <- "684380-99999"

weather_raw <- get_GSOD(
  station = kim_station,
  years   = 1993:2024
)
WeatherData <- weather_raw |>
  transmute(
    YEAR,
    MONTH,
    DAY,
    TEMP,
    MAX,
    MIN,
    PRCP,
    WDSP,
    MXSPD,
    I_HAIL
  ) |>
  arrange(YEAR, MONTH, DAY)

CHIRPS Rainfall for 2004

# coordinates  (24.81 E, 28.62 S)

lonlat <- data.frame(
lon = 24.81,
lat = -28.62
)

kimberley_chirps_2004 <- get_chirps(
object   = lonlat,
dates    = c("2004-01-01", "2004-12-31"),
server   = "ClimateSERV",
as.matrix = FALSE
)

Fetching data from ClimateSERV 

Getting your request...
glimpse(kimberley_chirps_2004)
Rows: 366
Columns: 5
$ id     <int> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, …
$ lon    <dbl> 24.81, 24.81, 24.81, 24.81, 24.81, 24.81, 24.81, 24.81, 24.81, …
$ lat    <dbl> -28.62, -28.62, -28.62, -28.62, -28.62, -28.62, -28.62, -28.62,…
$ date   <date> 2004-01-01, 2004-01-02, 2004-01-03, 2004-01-04, 2004-01-05, 20…
$ chirps <dbl> 0.000000, 0.000000, 0.000000, 0.000000, 12.762877, 21.295383, 1…

CHIRPS Rainfall for 2004 (Dronfield / Kimberley)

daily_chirps_2004 <- tibble(
date    = as.Date(kimberley_chirps_2004$date),
rain_mm = as.numeric(kimberley_chirps_2004$chirps)
)

head(daily_chirps_2004) 
# A tibble: 6 × 2
  date       rain_mm
  <date>       <dbl>
1 2004-01-01     0  
2 2004-01-02     0  
3 2004-01-03     0  
4 2004-01-04     0  
5 2004-01-05    12.8
6 2004-01-06    21.3
summary(daily_chirps_2004$rain_mm)
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
   0.00    0.00    0.00    1.26    0.00   28.76 
annual_chirps_2004 <- daily_chirps_2004 |>
mutate(YEAR = year(date)) |>
group_by(YEAR) |>
summarise(
total_rain_2004 = sum(rain_mm, na.rm = TRUE),
max_rain_2004   = max(rain_mm, na.rm = TRUE),
rain_days_2004  = sum(rain_mm > 0, na.rm = TRUE),
.groups = "drop"
)

annual_chirps_2004
# A tibble: 1 × 4
   YEAR total_rain_2004 max_rain_2004 rain_days_2004
  <dbl>           <dbl>         <dbl>          <int>
1  2004            461.          28.8             54
total_rain_2004 <- annual_chirps_2004$total_rain_2004
max_rain_2004   <- annual_chirps_2004$max_rain_2004
rain_days_2004  <- annual_chirps_2004$rain_days_2004

total_rain_2004
[1] 461.1258
max_rain_2004
[1] 28.75726
rain_days_2004
[1] 54

Weather Data

WeatherData_clean <- WeatherData |>
  filter(
    YEAR >= 1993,
    YEAR <= 2024
  ) |>
  mutate(
    date = make_date(YEAR, MONTH, DAY)
  ) |>
  # add CHIRPS daily rainfall into 2004 
  left_join(daily_chirps_2004, by = "date") |>
  mutate(
    PRCP = if_else(YEAR == 2004 & !is.na(rain_mm), rain_mm, PRCP)
  ) |>
  dplyr::select(-rain_mm) 
vultures.clean <- vultures |>
  clean_names() |>                               # ringing_date, laying_date
  rename_with(~ gsub("_", ".", .x)) |>           # ringing.date, laying.dater
  mutate(
    ringing.date = ringing.date |> as.character() |> str_trim(),
    ringing.date = if_else(ringing.date %in% c("", "NA", "?"),
                           NA_character_, ringing.date),
    ringing.date = ymd(gsub("/", "-", ringing.date)),
    
    laying.date  = laying.date |> as.character() |> str_trim(),
    laying.date  = if_else(laying.date %in% c("", "NA", "?"),
                           NA_character_, laying.date),
    laying.date  = ymd(gsub("/", "-", laying.date))
  )
Warning: There was 1 warning in `mutate()`.
ℹ In argument: `ringing.date = ymd(gsub("/", "-", ringing.date))`.
Caused by warning:
!  32 failed to parse.
glimpse(vultures.clean)
Rows: 2,657
Columns: 14
$ ringing.date            <date> 1993-09-21, 1993-09-21, 1993-09-19, 1993-09-1…
$ tree.drn                <chr> "1", "2", "3", "4", "5", "6", "8", "12", "13",…
$ southings               <chr> " 28 38 43.8", "28 38 48.2", "28 38 55.6", "28…
$ eastings                <chr> "24 47 14.5", "24 47 05.1", "24 47 15.4", "24 …
$ code                    <int> 1, 1, 1, 1, 1, 1, 1, 15, 1, 1, 1, 1, 3, 1, 5, …
$ ring                    <chr> "G19791", "G19792", "G19781", "G19780", "G1978…
$ colour.rings.r.top.down <chr> "M", "M", "M", "M", "M", "M", "M", "", "M", "M…
$ colour.rings.l.top.down <chr> "B, Y, Y", "B, Y, R", "B, W, G", "B, B, R", "B…
$ mass.g                  <int> 5900, 5000, 3900, 4950, 4950, 3800, 3550, NA, …
$ wing.length.mm          <int> 475, 425, 260, 310, 290, 275, 210, NA, 245, 29…
$ age.at.ringing.days     <int> 89, 80, 57, 64, 61, 59, 49, NA, 55, 61, 5, 92,…
$ value.at.ringing        <int> 34233, 34233, 34231, 34231, 34231, 34231, 3423…
$ hatching.date           <chr> "1993/06/24", "1993/07/03", "1993/07/24", "199…
$ laying.date             <date> 1993-04-29, 1993-05-08, 1993-05-29, 1993-05-2…
summary(vultures.clean)
  ringing.date          tree.drn          southings           eastings        
 Min.   :1993-09-19   Length:2657        Length:2657        Length:2657       
 1st Qu.:2004-06-06   Class :character   Class :character   Class :character  
 Median :2012-10-13   Mode  :character   Mode  :character   Mode  :character  
 Mean   :2011-08-31                                                           
 3rd Qu.:2019-08-21                                                           
 Max.   :2024-10-18                                                           
 NA's   :32                                                                   
      code             ring           colour.rings.r.top.down
 Min.   :  1.000   Length:2657        Length:2657            
 1st Qu.:  1.000   Class :character   Class :character       
 Median :  1.000   Mode  :character   Mode  :character       
 Mean   :  5.001                                             
 3rd Qu.:  6.000                                             
 Max.   :138.000                                             
 NA's   :4                                                   
 colour.rings.l.top.down     mass.g     wing.length.mm  age.at.ringing.days
 Length:2657             Min.   :   2   Min.   : 25.0   Min.   :  5.00     
 Class :character        1st Qu.:4200   1st Qu.:320.0   1st Qu.: 65.00     
 Mode  :character        Median :4800   Median :395.0   Median : 76.00     
                         Mean   :4610   Mean   :375.7   Mean   : 74.18     
                         3rd Qu.:5200   3rd Qu.:460.0   3rd Qu.: 86.00     
                         Max.   :6700   Max.   :570.0   Max.   :300.00     
                         NA's   :1225   NA's   :1192    NA's   :1198       
 value.at.ringing hatching.date       laying.date        
 Min.   :34231    Length:2657        Min.   :1900-01-02  
 1st Qu.:37541    Class :character   1st Qu.:2002-06-02  
 Median :40824    Mode  :character   Median :2011-05-13  
 Mean   :40489                       Mean   :2010-05-21  
 3rd Qu.:43386                       3rd Qu.:2018-05-26  
 Max.   :45583                       Max.   :2024-08-11  
 NA's   :1231                        NA's   :1197        
colSums(is.na(vultures.clean))
           ringing.date                tree.drn               southings 
                     32                       0                       0 
               eastings                    code                    ring 
                      0                       4                       0 
colour.rings.r.top.down colour.rings.l.top.down                  mass.g 
                      0                       0                    1225 
         wing.length.mm     age.at.ringing.days        value.at.ringing 
                   1192                    1198                    1231 
          hatching.date             laying.date 
                      0                    1197 

Active Nests per year

active.nests.raw <- vultures.clean |>
mutate(year = year(ringing.date)) |>
group_by(year) |>
summarise(active_nests = n(), .groups = "drop")

ggplot(active.nests.raw, aes(x = year, y = active_nests)) +
geom_col(fill = "#89c5d6", alpha = 0.85) +
labs(x = "Year", y = "Number of active nests") +
theme_bw() +
theme(panel.grid = element_blank())
Warning: Removed 1 row containing missing values or values outside the scale range
(`geom_col()`).

Success proportion

failed.nests.raw <- vultures.clean |>
mutate(year = year(ringing.date)) |>
group_by(year) |>
summarise(
active_nests  = n(),
success_nests = sum(!is.na(laying.date)),
success_prop  = success_nests / active_nests,
.groups       = "drop"
)

ggplot(failed.nests.raw, aes(x = success_prop)) +
geom_histogram(bins = 15, fill = "#89c5d6") +
labs(x = "Breeding success proportion", y = "Count") +
theme_bw() +
theme(panel.grid = element_blank())

Annual Weather Summaries

weather_yearly <- WeatherData_clean |>
group_by(YEAR) |>
summarise(
mean_temp  = mean(TEMP, na.rm = TRUE),
max_temp   = max(MAX, na.rm = TRUE),
min_temp   = min(MIN, na.rm = TRUE),
total_rain = sum(PRCP, na.rm = TRUE),
max_rain   = max(PRCP, na.rm = TRUE),
rain_days  = sum(PRCP > 0, na.rm = TRUE),
mean_wind  = mean(WDSP, na.rm = TRUE),
max_wind   = max(MXSPD, na.rm = TRUE),
hail_days  = sum(I_HAIL, na.rm = TRUE),
.groups    = "drop"
)

# Check the corrected 2004 rainfall

weather_yearly[weather_yearly$YEAR == 2004,
c("YEAR", "total_rain", "max_rain", "rain_days")]
# A tibble: 1 × 4
   YEAR total_rain max_rain rain_days
  <int>      <dbl>    <dbl>     <int>
1  2004       461.     28.8        54

Plot Annual Rainfall

ggplot(weather_yearly, aes(x = YEAR, y = total_rain)) +
geom_line() +
geom_point() +
theme_bw() +
labs(x = "Year", y = "Total annual rainfall (mm)")

Rain over time

library(grid) 

rain_yearly <- weather_yearly

rain_long <- rain_yearly |>
pivot_longer(
cols = c(total_rain, max_rain, rain_days),
names_to = "metric",
values_to = "value"
) |>
mutate(
metric = factor(
metric,
levels = c("total_rain", "max_rain", "rain_days")
),
metric_label = case_when(
metric == "total_rain" ~ "Total annual rainfall (mm)",
metric == "max_rain"   ~ "Max daily rainfall (mm)",
metric == "rain_days"  ~ "Number of rainy days (> 0 mm)"
),
metric_label = factor(
metric_label,
levels = c(
"Total annual rainfall (mm)",
"Max daily rainfall (mm)",
"Number of rainy days (> 0 mm)"
)
)
)

ggplot(rain_long, aes(x = YEAR, y = value)) +
geom_line(linewidth = 0.7) +
geom_point(size = 1.3) +
facet_wrap(~ metric_label, ncol = 1, scales = "free_y") +
theme_classic() +
labs(
x = "Year",
y = NULL
) +
theme(
strip.text    = element_text(size = 9, face = "bold"),
axis.title.x  = element_text(size = 10),
axis.text     = element_text(size = 8),
panel.spacing = unit(0.7, "lines")
)

Check min temp

# Yearly 

ggplot(weather_yearly, aes(x = YEAR, y = min_temp)) +
geom_line() +
geom_point() +
theme_bw() +
labs(x = "Year", y = "Annual minimum temperature (°C)")

# Daily 

weather_daily <- WeatherData_clean |>
group_by(YEAR, MONTH, DAY) |>
summarise(
mean_temp  = mean(TEMP, na.rm = TRUE),
max_temp   = max(MAX,  na.rm = TRUE),
min_temp   = min(MIN,  na.rm = TRUE),
total_rain = sum(PRCP, na.rm = TRUE),
max_rain   = max(PRCP, na.rm = TRUE),
rain_days  = sum(PRCP > 0, na.rm = TRUE),
mean_wind  = mean(WDSP, na.rm = TRUE),
max_wind   = max(MXSPD, na.rm = TRUE),
hail_days  = sum(I_HAIL, na.rm = TRUE),
.groups    = "drop"
)
Warning: There were 291 warnings in `summarise()`.
The first warning was:
ℹ In argument: `max_temp = max(MAX, na.rm = TRUE)`.
ℹ In group 9935: `YEAR = 2020`, `MONTH = 6`, `DAY = 15`.
Caused by warning in `max()`:
! no non-missing arguments to max; returning -Inf
ℹ Run `dplyr::last_dplyr_warnings()` to see the 290 remaining warnings.
temp_daily <- weather_daily |>
mutate(
min_temp = ifelse(is.infinite(min_temp), NA_real_, min_temp),
max_temp = ifelse(is.infinite(max_temp), NA_real_, max_temp)
)

summary(temp_daily$min_temp)
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
 -9.900   4.500  10.750   9.818  15.500  30.000       9 
summary(temp_daily$max_temp)
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
   5.50   23.70   28.60   28.18   32.90   45.70       1 

Yearly temp summaries

weather_yearly_temp <- weather_daily |>
filter(YEAR >= 1993, YEAR <= 2024) |>
group_by(YEAR) |>
summarise(
mean_temp = mean(mean_temp, na.rm = TRUE),
max_temp  = max(max_temp, na.rm = TRUE),
min_temp  = min(min_temp, na.rm = TRUE),
.groups   = "drop"
)

summary(weather_yearly_temp$min_temp)
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
 -9.900  -7.425  -6.450  -6.356  -5.450  -3.300 
summary(weather_yearly_temp$max_temp)
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  37.20   38.35   39.50   40.01   41.17   45.70 

Temps over time

min_temp <- ggplot(weather_yearly_temp, aes(YEAR, min_temp)) +
geom_line() +
geom_point() +
theme_bw() +
labs(x = "Year", y = "Annual minimum temperature (°C)")

max_temp <- ggplot(weather_yearly_temp, aes(x = YEAR, y = max_temp)) +
geom_line() +
geom_point() +
theme_bw() +
labs(
x = "Year",
y = "Annual maximum temperature (°C)",
title = "Annual Maximum Temperature (1993–2024)"
)

mean_temp <- ggplot(weather_yearly_temp, aes(x = YEAR, y = mean_temp)) +
geom_line() +
geom_point() +
theme_bw() +
labs(
x = "Year",
y = "Annual mean temperature (°C)",
title = "Annual Mean Temperature (1993–2024)"
)
temp_yearly_long <- weather_yearly_temp |>
pivot_longer(
cols = c(max_temp, min_temp, mean_temp),
names_to = "type",
values_to = "value"
) |>
mutate(
type = factor(
type,
levels = c("max_temp", "min_temp", "mean_temp"),
labels = c("Annual Maximum Temperature (°C)",
"Annual Minimum Temperature (°C)",
"Annual Mean Temperature (°C)")
)
)

ggplot(temp_yearly_long, aes(x = YEAR, y = value)) +
geom_line(size = 0.7) +
geom_point(size = 1.5) +
facet_wrap(~ type, scales = "free_y", ncol = 1) +
theme_classic(base_size = 11) +
theme(
strip.text   = element_text(size = 10, face = "bold"),
axis.title   = element_text(size = 9),
axis.text    = element_text(size = 8),
panel.spacing = unit(0.8, "lines")
) +
labs(
x = "Year",
y = "",
title = "Annual Temperature Trends (1993–2024)"
)
Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
ℹ Please use `linewidth` instead.

Lagged Weather summaries

weather_yearly_lag <- weather_yearly |>
arrange(YEAR) |>
mutate(
across(
c(mean_temp, max_temp, min_temp, total_rain, max_rain,
rain_days, mean_wind, max_wind, hail_days),
~ lag(.x, 1), .names = "lag_{.col}"
)
)

Daily weather with lagged variables

weather_lag <- WeatherData_clean |>
mutate(
mean_temp = readr::parse_number(as.character(TEMP)),
max_temp  = readr::parse_number(as.character(MAX)),
min_temp  = readr::parse_number(as.character(MIN)),
mean_wind = readr::parse_number(as.character(WDSP)),
max_wind  = readr::parse_number(as.character(MXSPD)),
rain      = readr::parse_number(as.character(PRCP)),
rain_day  = as.integer(replace_na(rain, 0) > 0),
hail      = I_HAIL
) |>
arrange(date)

Check for weather correlation

weather_yearly |>
  dplyr::select(
    mean_temp, max_temp, min_temp,
    total_rain, max_rain, rain_days,
    mean_wind, max_wind, hail_days
  )
# A tibble: 32 × 9
   mean_temp max_temp min_temp total_rain max_rain rain_days mean_wind max_wind
       <dbl>    <dbl>    <dbl>      <dbl>    <dbl>     <int>     <dbl>    <dbl>
 1      18.4     39.1     -7.8       568.    111          67      4.05     15.4
 2      17.8     39.6     -8         378.     37.3        39      3.75     13.9
 3      18.5     39.1     -7.2       567.     81.3        70      3.79     18  
 4      17.7     37.5     -6.5       516.     39.1        79      3.98     14.9
 5      18.5     40.9     -3.3       381.     33.0        65      3.76     13.9
 6      18.7     37.8     -6         550.     32          84      4.00     20.6
 7      19.4     37.9     -5         354.     42.2        62      3.95     28.3
 8      18.0     37.2     -8.1       578.     78.5        88      3.61     20.6
 9      19.0     39.9     -6.4       548.     27.4       105      3.78     27.8
10      19.4     43.6     -6.6       579.     44.7        83      3.71     24.2
# ℹ 22 more rows
# ℹ 1 more variable: hail_days <dbl>
cor_mat <- weather_yearly |>
dplyr::select(mean_temp, max_temp, min_temp, total_rain, max_rain,
rain_days, mean_wind, max_wind, hail_days) |>
cor(use = "pairwise.complete.obs")

cor_df <- as.data.frame(cor_mat) |>
rownames_to_column("var1") |>
pivot_longer(-var1, names_to = "var2", values_to = "cor")

ggplot(cor_df, aes(x = var1, y = var2, fill = cor)) +
geom_tile() +
scale_fill_gradient2(
low = "#b2182b", mid = "white", high = "#2166ac",
midpoint = 0, limits = c(-1, 1)
) +
labs(x = "", y = "", fill = "Correlation") +
theme_bw() +
theme(
axis.text.x = element_text(angle = 45, hjust = 1),
panel.grid  = element_blank()
)

All weather across years

weather_yearly |>
pivot_longer(-YEAR) |>
ggplot(aes(x = YEAR, y = value)) +
geom_line(colour = "#89c5d6") +
facet_wrap(~ name, scales = "free_y") +
labs(x = "Year", y = "Value") +
theme_bw() +
theme(panel.grid = element_blank())

# Effort Data 

active.nests <- vultures.clean |>
mutate(year = year(ringing.date)) |>
group_by(year) |>
summarise(active_nests = n(), .groups = "drop")

effort_weather <- active.nests |>
left_join(weather_yearly, by = c("year" = "YEAR"))

Year only model

m_year <- glm.nb(active_nests ~ year, data = active.nests)
summary(m_year)

Call:
glm.nb(formula = active_nests ~ year, data = active.nests, init.theta = 283.4337529, 
    link = log)

Coefficients:
              Estimate Std. Error z value Pr(>|z|)    
(Intercept) -53.188260   4.927882  -10.79   <2e-16 ***
year          0.028659   0.002451   11.69   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for Negative Binomial(283.4338) family taken to be 1)

    Null deviance: 173.368  on 31  degrees of freedom
Residual deviance:  34.003  on 30  degrees of freedom
  (1 observation deleted due to missingness)
AIC: 246.56

Number of Fisher Scoring iterations: 1

              Theta:  283 
          Std. Err.:  329 

 2 x log-likelihood:  -240.558 
check_overdispersion(m_year)
[1] 1.067024

Autocorrelation in residuals

res_effort <- residuals(m_year, type = "pearson")

acf(res_effort,
main = "ACF of residuals: breeding effort (active nests)")

Rainfall

effort_weather <- effort_weather |>
mutate(
z_total_rain = scale(total_rain),
z_max_rain   = scale(max_rain),
z_rain_days  = scale(rain_days)
)

m_rain_only <- glm.nb(
active_nests ~ z_total_rain + z_max_rain + z_rain_days,
data = effort_weather
)
summary(m_rain_only)

Call:
glm.nb(formula = active_nests ~ z_total_rain + z_max_rain + z_rain_days, 
    data = effort_weather, init.theta = 18.0307885, link = log)

Coefficients:
             Estimate Std. Error z value Pr(>|z|)    
(Intercept)   4.39736    0.04606  95.461   <2e-16 ***
z_total_rain -0.18087    0.07301  -2.477   0.0132 *  
z_max_rain   -0.01308    0.05642  -0.232   0.8167    
z_rain_days   0.11469    0.06637   1.728   0.0840 .  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for Negative Binomial(18.0308) family taken to be 1)

    Null deviance: 41.511  on 31  degrees of freedom
Residual deviance: 32.342  on 28  degrees of freedom
  (1 observation deleted due to missingness)
AIC: 294.97

Number of Fisher Scoring iterations: 1

              Theta:  18.03 
          Std. Err.:  5.51 

 2 x log-likelihood:  -284.973 

Temp

effort_weather <- effort_weather |>
mutate(
z_mean_temp = scale(mean_temp),
z_max_temp  = scale(max_temp),
z_min_temp  = scale(min_temp)
)

max_temp_mod <- glm.nb(
active_nests ~ z_max_temp,
data = effort_weather
)
summary(max_temp_mod)

Call:
glm.nb(formula = active_nests ~ z_max_temp, data = effort_weather, 
    init.theta = 15.77748683, link = log)

Coefficients:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept)  4.40154    0.04864  90.494   <2e-16 ***
z_max_temp   0.10583    0.04911   2.155   0.0312 *  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for Negative Binomial(15.7775) family taken to be 1)

    Null deviance: 37.216  on 31  degrees of freedom
Residual deviance: 32.494  on 30  degrees of freedom
  (1 observation deleted due to missingness)
AIC: 294.66

Number of Fisher Scoring iterations: 1

              Theta:  15.78 
          Std. Err.:  4.72 

 2 x log-likelihood:  -288.663 
min_temp_mod <- glm.nb(
active_nests ~ z_min_temp,
data = effort_weather
)
summary(min_temp_mod)

Call:
glm.nb(formula = active_nests ~ z_min_temp, data = effort_weather, 
    init.theta = 13.63264473, link = log)

Coefficients:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept)  4.40648    0.05171  85.218   <2e-16 ***
z_min_temp   0.03703    0.05254   0.705    0.481    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for Negative Binomial(13.6326) family taken to be 1)

    Null deviance: 32.929  on 31  degrees of freedom
Residual deviance: 32.466  on 30  degrees of freedom
  (1 observation deleted due to missingness)
AIC: 298.61

Number of Fisher Scoring iterations: 1

              Theta:  13.63 
          Std. Err.:  3.97 

 2 x log-likelihood:  -292.606 
mean_temp_mod <- glm.nb(
active_nests ~ z_mean_temp,
data = effort_weather
)
summary(mean_temp_mod)

Call:
glm.nb(formula = active_nests ~ z_mean_temp, data = effort_weather, 
    init.theta = 13.63790192, link = log)

Coefficients:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept)  4.40648    0.05170  85.233   <2e-16 ***
z_mean_temp  0.03478    0.05251   0.662    0.508    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for Negative Binomial(13.6379) family taken to be 1)

    Null deviance: 32.939  on 31  degrees of freedom
Residual deviance: 32.476  on 30  degrees of freedom
  (1 observation deleted due to missingness)
AIC: 298.61

Number of Fisher Scoring iterations: 1

              Theta:  13.64 
          Std. Err.:  3.97 

 2 x log-likelihood:  -292.605 

Wind

effort_weather <- effort_weather |>
mutate(
z_mean_wind = scale(mean_wind),
z_max_wind  = scale(max_wind)
)

m_wind_only <- glm.nb(
active_nests ~ z_mean_wind + z_max_wind,
data = effort_weather
)
summary(m_wind_only)

Call:
glm.nb(formula = active_nests ~ z_mean_wind + z_max_wind, data = effort_weather, 
    init.theta = 15.33069655, link = log)

Coefficients:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept)  4.40243    0.04922  89.438   <2e-16 ***
z_mean_wind -0.06706    0.05215  -1.286    0.199    
z_max_wind  -0.05428    0.05231  -1.038    0.299    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for Negative Binomial(15.3307) family taken to be 1)

    Null deviance: 36.340  on 31  degrees of freedom
Residual deviance: 32.489  on 29  degrees of freedom
  (1 observation deleted due to missingness)
AIC: 297.43

Number of Fisher Scoring iterations: 1

              Theta:  15.33 
          Std. Err.:  4.56 

 2 x log-likelihood:  -289.431 

Hail

effort_weather <- effort_weather |>
mutate(z_hail_days = scale(hail_days))

m_hail_only <- glm.nb(
active_nests ~ z_hail_days,
data = effort_weather
)
summary(m_hail_only)

Call:
glm.nb(formula = active_nests ~ z_hail_days, data = effort_weather, 
    init.theta = 15.4472834, link = log)

Coefficients:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept)  4.40226    0.04907  89.721   <2e-16 ***
z_hail_days  0.09695    0.04926   1.968    0.049 *  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for Negative Binomial(15.4473) family taken to be 1)

    Null deviance: 36.569  on 31  degrees of freedom
Residual deviance: 32.511  on 30  degrees of freedom
  (1 observation deleted due to missingness)
AIC: 295.25

Number of Fisher Scoring iterations: 1

              Theta:  15.45 
          Std. Err.:  4.60 

 2 x log-likelihood:  -289.249 

Active nests across years

active.nests <- vultures.clean |>
mutate(year = year(ringing.date)) |>
group_by(year) |>
summarise(active_nests = n(), .groups = "drop")

m_year <- glm.nb(active_nests ~ year, data = active.nests)

active.nests_plot <- active.nests |>
mutate(pred_year = predict(m_year, newdata = active.nests, type = "response"))

ggplot(active.nests_plot, aes(x = year, y = active_nests)) +
geom_col(fill = "#89c5d6", alpha = 0.85) +
geom_line(aes(y = pred_year), linewidth = 1.1, colour = "black") +
labs(
x = "Year",
y = "Number of active nests"
) +
theme_bw() +
theme(panel.grid = element_blank())
Warning: Removed 1 row containing missing values or values outside the scale range
(`geom_col()`).
Warning: Removed 1 row containing missing values or values outside the scale range
(`geom_line()`).

Effects of max temp on effort

effort_weather <- active.nests |>
left_join(weather_yearly, by = c("year" = "YEAR")) |>
mutate(
z_total_rain = as.numeric(scale(total_rain)),
z_max_rain   = as.numeric(scale(max_rain)),
z_rain_days  = as.numeric(scale(rain_days)),
z_mean_temp  = as.numeric(scale(mean_temp)),
z_max_temp   = as.numeric(scale(max_temp)),
z_min_temp   = as.numeric(scale(min_temp)),
z_mean_wind  = as.numeric(scale(mean_wind)),
z_max_wind   = as.numeric(scale(max_wind)),
z_hail_days  = as.numeric(scale(hail_days))
)

m_eff_maxtemp <- glm.nb(
active_nests ~ z_max_temp,
data = effort_weather
)
summary(m_eff_maxtemp)

Call:
glm.nb(formula = active_nests ~ z_max_temp, data = effort_weather, 
    init.theta = 15.77748683, link = log)

Coefficients:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept)  4.40154    0.04864  90.494   <2e-16 ***
z_max_temp   0.10583    0.04911   2.155   0.0312 *  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for Negative Binomial(15.7775) family taken to be 1)

    Null deviance: 37.216  on 31  degrees of freedom
Residual deviance: 32.494  on 30  degrees of freedom
  (1 observation deleted due to missingness)
AIC: 294.66

Number of Fisher Scoring iterations: 1

              Theta:  15.78 
          Std. Err.:  4.72 

 2 x log-likelihood:  -288.663 
new_eff_temp <- tibble(
z_max_temp = seq(
min(effort_weather$z_max_temp, na.rm = TRUE),
max(effort_weather$z_max_temp, na.rm = TRUE),
length.out = 100
)
) |>
mutate(
pred_nests = predict(
m_eff_maxtemp,
newdata = cur_data(),
type = "response"
)
)
Warning: There was 1 warning in `mutate()`.
ℹ In argument: `pred_nests = predict(m_eff_maxtemp, newdata = cur_data(), type
  = "response")`.
Caused by warning:
! `cur_data()` was deprecated in dplyr 1.1.0.
ℹ Please use `pick()` instead.
ggplot(effort_weather, aes(x = z_max_temp, y = active_nests)) +
geom_point(alpha = 0.7) +
geom_line(
data = new_eff_temp,
aes(x = z_max_temp, y = pred_nests),
linewidth = 1.1
) +
labs(
x = "Annual maximum temperature (scaled)",
y = "Number of active nests"
) +
theme_bw() +
theme(panel.grid = element_blank())
Warning: Removed 1 row containing missing values or values outside the scale range
(`geom_point()`).

acf(active.nests$active_nests)

Success DF

failed.nests <- vultures.clean |>
mutate(year = year(ringing.date)) |>
filter(!is.na(year)) |>
group_by(year) |>
summarise(
active_nests  = n(),
success_nests = sum(!is.na(laying.date)),
failed_nests  = active_nests - success_nests,
success_prop  = success_nests / active_nests,
.groups       = "drop"
) |>
filter(active_nests > 0)

Year only model

success_year <- glm(
cbind(success_nests, failed_nests) ~ year,
family = binomial,
data   = failed.nests
)
summary(success_year)

Call:
glm(formula = cbind(success_nests, failed_nests) ~ year, family = binomial, 
    data = failed.nests)

Coefficients:
             Estimate Std. Error z value Pr(>|z|)    
(Intercept) 45.434641   8.858715   5.129 2.92e-07 ***
year        -0.022505   0.004405  -5.109 3.24e-07 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 107.689  on 31  degrees of freedom
Residual deviance:  81.283  on 30  degrees of freedom
AIC: 237.58

Number of Fisher Scoring iterations: 3
check_overdispersion(success_year)
[1] 2.670817

Rain

success_df <- failed.nests |>
left_join(weather_yearly, by = c("year" = "YEAR")) |>
mutate(
z_total_rain = as.numeric(scale(total_rain)),
z_max_rain   = as.numeric(scale(max_rain)),
z_rain_days  = as.numeric(scale(rain_days))
)

success_rain <- glm(
cbind(success_nests, failed_nests) ~
z_total_rain + z_max_rain + z_rain_days,
family = binomial,
data   = success_df
)
summary(success_rain)

Call:
glm(formula = cbind(success_nests, failed_nests) ~ z_total_rain + 
    z_max_rain + z_rain_days, family = binomial, data = success_df)

Coefficients:
             Estimate Std. Error z value Pr(>|z|)    
(Intercept)   0.18710    0.03966   4.717 2.39e-06 ***
z_total_rain  0.08304    0.06426   1.292   0.1963    
z_max_rain    0.02901    0.05086   0.570   0.5685    
z_rain_days  -0.09864    0.05967  -1.653   0.0983 .  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 107.69  on 31  degrees of freedom
Residual deviance: 102.75  on 28  degrees of freedom
AIC: 263.05

Number of Fisher Scoring iterations: 3
check_overdispersion(success_rain)
[1] 3.581363

Temp

success_df <- success_df |>
mutate(
z_mean_temp = as.numeric(scale(mean_temp)),
z_max_temp  = as.numeric(scale(max_temp)),
z_min_temp  = as.numeric(scale(min_temp))
)

success_minT <- glm(
cbind(success_nests, failed_nests) ~ z_min_temp,
family = binomial,
data   = success_df
)
summary(success_minT)

Call:
glm(formula = cbind(success_nests, failed_nests) ~ z_min_temp, 
    family = binomial, data = success_df)

Coefficients:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept)  0.17203    0.03931   4.376 1.21e-05 ***
z_min_temp   0.16071    0.04191   3.834 0.000126 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 107.689  on 31  degrees of freedom
Residual deviance:  92.877  on 30  degrees of freedom
AIC: 249.18

Number of Fisher Scoring iterations: 3
check_overdispersion(success_minT)
[1] 3.031976
success_maxT <- glm(
cbind(success_nests, failed_nests) ~ z_max_temp,
family = binomial,
data   = success_df
)
summary(success_maxT)

Call:
glm(formula = cbind(success_nests, failed_nests) ~ z_max_temp, 
    family = binomial, data = success_df)

Coefficients:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept)  0.18310    0.03944   4.643 3.44e-06 ***
z_max_temp  -0.06103    0.03839  -1.590    0.112    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 107.69  on 31  degrees of freedom
Residual deviance: 105.16  on 30  degrees of freedom
AIC: 261.46

Number of Fisher Scoring iterations: 3
check_overdispersion(success_maxT)
[1] 3.433179
success_meanT <- glm(
cbind(success_nests, failed_nests) ~ z_mean_temp,
family = binomial,
data   = success_df
)
summary(success_meanT)

Call:
glm(formula = cbind(success_nests, failed_nests) ~ z_mean_temp, 
    family = binomial, data = success_df)

Coefficients:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept)  0.17461    0.03922   4.452 8.51e-06 ***
z_mean_temp  0.05587    0.03856   1.449    0.147    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 107.69  on 31  degrees of freedom
Residual deviance: 105.59  on 30  degrees of freedom
AIC: 261.89

Number of Fisher Scoring iterations: 3
check_overdispersion(success_meanT)
[1] 3.4372

WInd

success_df <- success_df |>
mutate(
z_mean_wind = as.numeric(scale(mean_wind)),
z_max_wind  = as.numeric(scale(max_wind))
)

m_success_wind <- glm(
cbind(success_nests, failed_nests) ~ z_mean_wind + z_max_wind,
family = binomial,
data   = success_df
)
summary(m_success_wind)

Call:
glm(formula = cbind(success_nests, failed_nests) ~ z_mean_wind + 
    z_max_wind, family = binomial, data = success_df)

Coefficients:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept)  0.17646    0.03937   4.482  7.4e-06 ***
z_mean_wind  0.01436    0.04176   0.344    0.731    
z_max_wind  -0.01688    0.04351  -0.388    0.698    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 107.69  on 31  degrees of freedom
Residual deviance: 107.48  on 29  degrees of freedom
AIC: 265.78

Number of Fisher Scoring iterations: 3
check_overdispersion(m_success_wind)
[1] 3.625179

Hail

success_df <- success_df |>
mutate(
z_hail_days = as.numeric(scale(hail_days))
)

m_success_hail <- glm(
cbind(success_nests, failed_nests) ~ z_hail_days,
family = binomial,
data   = success_df
)
summary(m_success_hail)

Call:
glm(formula = cbind(success_nests, failed_nests) ~ z_hail_days, 
    family = binomial, data = success_df)

Coefficients:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept)  0.18275    0.03940   4.638 3.52e-06 ***
z_hail_days -0.05984    0.03647  -1.641    0.101    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 107.69  on 31  degrees of freedom
Residual deviance: 105.00  on 30  degrees of freedom
AIC: 261.3

Number of Fisher Scoring iterations: 3
check_overdispersion(m_success_hail)
[1] 3.424451

60 and 30 day pre laying period

timing_df <- vultures.clean |>
filter(!is.na(laying.date)) |>
mutate(
year    = year(laying.date),
lay_DOY = yday(laying.date)
)

summarise_60 <- function(lay_date, weather) {
w <- weather |>
filter(date > (lay_date - days(60)), date <= lay_date)

tibble(
mean_temp_60  = mean(w$mean_temp, na.rm = TRUE),
max_temp_60   = mean(w$max_temp, na.rm = TRUE),
min_temp_60   = mean(w$min_temp, na.rm = TRUE),
total_rain_60 = sum(w$rain, na.rm = TRUE),
rain_days_60  = sum(w$rain_day, na.rm = TRUE),
mean_wind_60  = mean(w$mean_wind, na.rm = TRUE),
max_wind_60   = mean(w$max_wind, na.rm = TRUE),
hail_days_60  = sum(w$hail, na.rm = TRUE)
)
}

timing_weather_60 <- timing_df |>
mutate(win60 = purrr::map(laying.date, summarise_60, weather = weather_lag)) |>
unnest(win60)

summarise_30 <- function(lay_date, weather) {
  w <- weather |>
    filter(date > (lay_date - days(30)), date <= lay_date)

  tibble(
    mean_temp_30  = mean(w$mean_temp, na.rm = TRUE),
    max_temp_30   = mean(w$max_temp,  na.rm = TRUE),
    min_temp_30   = mean(w$min_temp,  na.rm = TRUE),
    total_rain_30 = sum(w$rain,       na.rm = TRUE),
    rain_days_30  = sum(w$rain_day,   na.rm = TRUE),
    mean_wind_30  = mean(w$mean_wind, na.rm = TRUE),
    max_wind_30   = mean(w$max_wind,  na.rm = TRUE),
    hail_days_30  = sum(w$hail,       na.rm = TRUE)
  )
}

timing_weather_30 <- timing_df |>
  mutate(win30 = purrr::map(laying.date, summarise_30, weather = weather_lag)) |>
  tidyr::unnest(win30)

Temp

m_timing_temp_60 <- lm(
lay_DOY ~ mean_temp_60 + max_temp_60 + min_temp_60,
data = timing_weather_60
)
summary(m_timing_temp_60)

Call:
lm(formula = lay_DOY ~ mean_temp_60 + max_temp_60 + min_temp_60, 
    data = timing_weather_60)

Residuals:
    Min      1Q  Median      3Q     Max 
-88.511  -5.887  -1.339   4.323 145.827 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept)  230.2640     4.2240  54.514  < 2e-16 ***
mean_temp_60   2.6343     0.6348   4.150 3.52e-05 ***
max_temp_60   -3.0285     0.3587  -8.443  < 2e-16 ***
min_temp_60   -5.5674     0.3681 -15.123  < 2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 10.13 on 1455 degrees of freedom
  (1 observation deleted due to missingness)
Multiple R-squared:  0.6752,    Adjusted R-squared:  0.6745 
F-statistic:  1008 on 3 and 1455 DF,  p-value: < 2.2e-16

Rain

m_timing_rain_60 <- lm(
lay_DOY ~ total_rain_60 + rain_days_60,
data = timing_weather_60
)
summary(m_timing_rain_60)

Call:
lm(formula = lay_DOY ~ total_rain_60 + rain_days_60, data = timing_weather_60)

Residuals:
     Min       1Q   Median       3Q      Max 
-166.729   -9.039   -2.407    7.930  134.078 

Coefficients:
               Estimate Std. Error t value Pr(>|t|)    
(Intercept)   168.72873    0.85555 197.217  < 2e-16 ***
total_rain_60  -0.06469    0.01242  -5.207 2.19e-07 ***
rain_days_60   -0.96694    0.09618 -10.054  < 2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 15.69 on 1457 degrees of freedom
Multiple R-squared:  0.2568,    Adjusted R-squared:  0.2557 
F-statistic: 251.7 on 2 and 1457 DF,  p-value: < 2.2e-16

Hail

m_timing_rain_60 <- lm(
lay_DOY ~ total_rain_60 + rain_days_60,
data = timing_weather_60
)
summary(m_timing_rain_60)

Call:
lm(formula = lay_DOY ~ total_rain_60 + rain_days_60, data = timing_weather_60)

Residuals:
     Min       1Q   Median       3Q      Max 
-166.729   -9.039   -2.407    7.930  134.078 

Coefficients:
               Estimate Std. Error t value Pr(>|t|)    
(Intercept)   168.72873    0.85555 197.217  < 2e-16 ***
total_rain_60  -0.06469    0.01242  -5.207 2.19e-07 ***
rain_days_60   -0.96694    0.09618 -10.054  < 2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 15.69 on 1457 degrees of freedom
Multiple R-squared:  0.2568,    Adjusted R-squared:  0.2557 
F-statistic: 251.7 on 2 and 1457 DF,  p-value: < 2.2e-16

Test

plot(failed.nests$year, failed.nests$success_prop)


Objective 1

objective1_table <- tribble(
  ~`Hypothesis / justification`, ~Model,

  # ---- Temperature ----
  "**Temp**", NA,
  "Warmer years increase energetic cost of adult birds, reducing likelihood to breed.",
  "active_nests ~ mean_temp",

  "Extreme heat events cause thermal stress, reducing breeding effort (active cooling is energetically costly and may reduce adult condition).",
  "active_nests ~ max_temp",

  "Cold extremes increase thermoregulatory demand (often increasing time spent foraging), reducing effort put into breeding",
  "active_nests ~ min_temp",

  # ---- Rain ----
  "**Rain**", NA,
  "High rainfall reduces foraging efficiency and adult condition, reducing breeding effort.",
  "active_nests ~ total_rain",

  "Prolonged rainfall conditions reduce breeding effort by limiting foraging opportunity (fewer or weaker thermals for soaring)",
  "active_nests ~ rain_days",

  "Intense rainfall events decrease foraging efficiency and nest attendance, reducing breeding effort.",
  "active_nests ~ max_rain",

  "Previous year’s rainfall influences food availability and adult condition, affecting breeding effort in the following year.",
  "active_nests ~ lag_total_rain",

  # ---- Wind ----
  "**Wind**", NA,
  "Persistent windy conditions increase the cost of flying, limiting foraging and the energy available for reproduction, reducing breeding effort",
  "active_nests ~ mean_wind",

  "Extreme winds increase nest damage risk, reducing the likelihood that adults will attempt breeding (WBV often reuse nests from previous years) ",
  "active_nests ~ max_wind",

  # ---- Hail ----
  "**Hail**", NA,
  "Severe storm events increase disturbance and nest damage risk, reducing breeding effort",
  "active_nests ~ hail_days",

  # ---- Joint Effects ----
  "**Joint effects**", NA,
  "Cold stress combined with prolonged wet conditions jointly limit breeding effort",
  "active_nests ~ min_temp + total_rain",

  "Cold stress combined with intense rainfall events reduces breeding effort",
  "active_nests ~ min_temp + max_rain",

  "Cold extremes combined with extreme wind increase energetic costs and breeding risk, reducing breeding effort",
  "active_nests ~ min_temp + max_wind",

  "Severe storm exposure combined with extreme wind increases nest disturbance, reducing breeding effort",
  "active_nests ~ hail_days + max_wind",

  "Reduced thermals for soaring combined with high flight costs limit breeding effort",
  "active_nests ~ max_rain + max_wind",
  
  "Heat stress combined with prolonged wet conditions jointly reduce breeding effort by increasing energetic costs while limiting foraging",
  "active_nests ~ max_temp + rain_days",
  
 " High annual rainfall and frequent hail events each impose energetic and disturbance costs that reduce breeding effort.",
"active_nests ~ hail_days + rain_days",

  # ---- Interactive Effects ----
  "**Interactive effects**", NA,
  "The negative effect of cold extremes on breeding effort is increased in persistently wet years.",
  "active_nests ~ min_temp * rain_days",

  "Cold-related energetic stress is amplified in windy years due to increased heat loss and flight costs",
  "active_nests ~ min_temp * mean_wind",

  "Adult condition resulting from the previous year interacts with current-year heat stress to influence breeding effort",
  "active_nests ~ lag_total_rain * max_temp"
)

objective1_table_clean <- objective1_table |>
  dplyr::mutate(Model = ifelse(is.na(Model), "", Model))

kable(
  objective1_table_clean,
  align = c("l", "l"),
  caption = "Objective 1 candidate model set for breeding effort (active nests)"
)
Objective 1 candidate model set for breeding effort (active nests)
Hypothesis / justification Model
Temp
Warmer years increase energetic cost of adult birds, reducing likelihood to breed. active_nests ~ mean_temp
Extreme heat events cause thermal stress, reducing breeding effort (active cooling is energetically costly and may reduce adult condition). active_nests ~ max_temp
Cold extremes increase thermoregulatory demand (often increasing time spent foraging), reducing effort put into breeding active_nests ~ min_temp
Rain
High rainfall reduces foraging efficiency and adult condition, reducing breeding effort. active_nests ~ total_rain
Prolonged rainfall conditions reduce breeding effort by limiting foraging opportunity (fewer or weaker thermals for soaring) active_nests ~ rain_days
Intense rainfall events decrease foraging efficiency and nest attendance, reducing breeding effort. active_nests ~ max_rain
Previous year’s rainfall influences food availability and adult condition, affecting breeding effort in the following year. active_nests ~ lag_total_rain
Wind
Persistent windy conditions increase the cost of flying, limiting foraging and the energy available for reproduction, reducing breeding effort active_nests ~ mean_wind
Extreme winds increase nest damage risk, reducing the likelihood that adults will attempt breeding (WBV often reuse nests from previous years) active_nests ~ max_wind
Hail
Severe storm events increase disturbance and nest damage risk, reducing breeding effort active_nests ~ hail_days
Joint effects
Cold stress combined with prolonged wet conditions jointly limit breeding effort active_nests ~ min_temp + total_rain
Cold stress combined with intense rainfall events reduces breeding effort active_nests ~ min_temp + max_rain
Cold extremes combined with extreme wind increase energetic costs and breeding risk, reducing breeding effort active_nests ~ min_temp + max_wind
Severe storm exposure combined with extreme wind increases nest disturbance, reducing breeding effort active_nests ~ hail_days + max_wind
Reduced thermals for soaring combined with high flight costs limit breeding effort active_nests ~ max_rain + max_wind
Heat stress combined with prolonged wet conditions jointly reduce breeding effort by increasing energetic costs while limiting foraging active_nests ~ max_temp + rain_days
High annual rainfall and frequent hail events each impose energetic and disturbance costs that reduce breeding effort. active_nests ~ hail_days + rain_days
Interactive effects
The negative effect of cold extremes on breeding effort is increased in persistently wet years. active_nests ~ min_temp * rain_days
Cold-related energetic stress is amplified in windy years due to increased heat loss and flight costs active_nests ~ min_temp * mean_wind
Adult condition resulting from the previous year interacts with current-year heat stress to influence breeding effort active_nests ~ lag_total_rain * max_temp
# Annual effort: number of active nests per year
active.nests <- vultures.clean |>
  mutate(year = year(ringing.date)) |>
  group_by(year) |>
  summarise(active_nests = n(), .groups = "drop")


# Join to weather, create log effort response
effort_weather <- active.nests |>
  left_join(weather_yearly_lag, by = c("year" = "YEAR")) |>
  filter(!is.na(lag_total_rain)) |>
  mutate(
    log_active_nests = log(active_nests)
  )
# checking residuals
par(mfrow = c(1, 2))

Fitingt then no-year models

# Baseline 
m_null <- lm(log_active_nests ~ 1, data = effort_weather)

# Single predictor models 
m_mean_temp      <- lm(log_active_nests ~ mean_temp,      data = effort_weather)
m_max_temp       <- lm(log_active_nests ~ max_temp,       data = effort_weather)
m_min_temp       <- lm(log_active_nests ~ min_temp,       data = effort_weather)

m_total_rain     <- lm(log_active_nests ~ total_rain,     data = effort_weather)
m_rain_days      <- lm(log_active_nests ~ rain_days,      data = effort_weather)
m_max_rain       <- lm(log_active_nests ~ max_rain,       data = effort_weather)
m_lag_total_rain <- lm(log_active_nests ~ lag_total_rain, data = effort_weather)

m_mean_wind      <- lm(log_active_nests ~ mean_wind,      data = effort_weather)
m_max_wind       <- lm(log_active_nests ~ max_wind,       data = effort_weather)
m_hail_days      <- lm(log_active_nests ~ hail_days,      data = effort_weather)

# Joint effects 
m_minT_totalR    <- lm(log_active_nests ~ min_temp + total_rain,  data = effort_weather)
m_minT_maxR      <- lm(log_active_nests ~ min_temp + max_rain,    data = effort_weather)
m_minT_maxW      <- lm(log_active_nests ~ min_temp + max_wind,    data = effort_weather)
m_hail_maxW      <- lm(log_active_nests ~ hail_days + max_wind,   data = effort_weather)
m_maxR_maxW      <- lm(log_active_nests ~ max_rain + max_wind,    data = effort_weather)
m_maxT_rainDays  <- lm(log_active_nests ~ max_temp + rain_days,   data = effort_weather)
m_hail_rainDays  <- lm(log_active_nests ~ hail_days + rain_days,  data = effort_weather)

# Interactions 
m_minT_x_rainDays <- lm(log_active_nests ~ min_temp * rain_days,      data = effort_weather)
m_minT_x_meanWind <- lm(log_active_nests ~ min_temp * mean_wind,      data = effort_weather)
m_lagR_x_maxT     <- lm(log_active_nests ~ lag_total_rain * max_temp, data = effort_weather)

Fit with year models

# Baseline 
m_year <- lm(log_active_nests ~ year, data = effort_weather)

# Single predictor models 
m_mean_temp_year      <- lm(log_active_nests ~ year + mean_temp,      data = effort_weather)
m_max_temp_year       <- lm(log_active_nests ~ year + max_temp,       data = effort_weather)
m_min_temp_year       <- lm(log_active_nests ~ year + min_temp,       data = effort_weather)

m_total_rain_year     <- lm(log_active_nests ~ year + total_rain,     data = effort_weather)
m_rain_days_year      <- lm(log_active_nests ~ year + rain_days,      data = effort_weather)
m_max_rain_year       <- lm(log_active_nests ~ year + max_rain,       data = effort_weather)
m_lag_total_rain_year <- lm(log_active_nests ~ year + lag_total_rain, data = effort_weather)

m_mean_wind_year      <- lm(log_active_nests ~ year + mean_wind,      data = effort_weather)
m_max_wind_year       <- lm(log_active_nests ~ year + max_wind,       data = effort_weather)
m_hail_days_year      <- lm(log_active_nests ~ year + hail_days,      data = effort_weather)

# Joint effects 
m_minT_totalR_year    <- lm(log_active_nests ~ year + min_temp + total_rain, data = effort_weather)
m_minT_maxR_year      <- lm(log_active_nests ~ year + min_temp + max_rain,   data = effort_weather)
m_minT_maxW_year      <- lm(log_active_nests ~ year + min_temp + max_wind,   data = effort_weather)
m_hail_maxW_year      <- lm(log_active_nests ~ year + hail_days + max_wind,  data = effort_weather)
m_maxR_maxW_year      <- lm(log_active_nests ~ year + max_rain + max_wind,   data = effort_weather)
m_maxT_rainDays_year  <- lm(log_active_nests ~ year + max_temp + rain_days,  data = effort_weather)
m_hail_rainDays_year  <- lm(log_active_nests ~ year + hail_days + rain_days, data = effort_weather)

# Interactions 
m_minT_x_rainDays_year <- lm(log_active_nests ~ year + min_temp * rain_days,      data = effort_weather)
m_minT_x_meanWind_year <- lm(log_active_nests ~ year + min_temp * mean_wind,      data = effort_weather)
m_lagR_x_maxT_year     <- lm(log_active_nests ~ year + lag_total_rain * max_temp, data = effort_weather)

Two separate AICc model selection runs

# AIC 1: without year
cand_effort_no_year <- list(
  null            = m_null,

  mean_temp       = m_mean_temp,
  max_temp        = m_max_temp,
  min_temp        = m_min_temp,
  total_rain      = m_total_rain,
  rain_days       = m_rain_days,
  max_rain        = m_max_rain,
  lag_total_rain  = m_lag_total_rain,
  mean_wind       = m_mean_wind,
  max_wind        = m_max_wind,
  hail_days       = m_hail_days,

  minT_totalR     = m_minT_totalR,
  minT_maxR       = m_minT_maxR,
  minT_maxW       = m_minT_maxW,
  hail_maxW       = m_hail_maxW,
  maxR_maxW       = m_maxR_maxW,
  maxT_rainDays   = m_maxT_rainDays,
  hail_rainDays   = m_hail_rainDays,

  minT_x_rainDays = m_minT_x_rainDays,
  minT_x_meanWind = m_minT_x_meanWind,
  lagR_x_maxT     = m_lagR_x_maxT
)

aic_effort_no_year <- aictab(cand.set = cand_effort_no_year, modnames = names(cand_effort_no_year))
aic_effort_no_year

Model selection based on AICc:

                K  AICc Delta_AICc AICcWt Cum.Wt    LL
hail_days       3 12.04       0.00   0.17   0.17 -2.57
max_temp        3 12.77       0.74   0.12   0.29 -2.94
total_rain      3 12.83       0.79   0.12   0.41 -2.97
hail_maxW       4 13.16       1.12   0.10   0.51 -1.81
null            2 14.26       2.22   0.06   0.57 -4.92
maxT_rainDays   4 14.37       2.33   0.05   0.62 -2.41
max_wind        3 14.39       2.35   0.05   0.68 -3.75
hail_rainDays   4 14.44       2.40   0.05   0.73 -2.45
mean_wind       3 14.84       2.81   0.04   0.77 -3.98
lag_total_rain  3 15.07       3.03   0.04   0.81 -4.09
minT_totalR     4 15.45       3.41   0.03   0.84 -2.96
lagR_x_maxT     5 15.71       3.67   0.03   0.87 -1.65
max_rain        3 15.74       3.70   0.03   0.89 -4.43
maxR_maxW       4 15.84       3.81   0.03   0.92 -3.15
mean_temp       3 16.46       4.42   0.02   0.94 -4.79
min_temp        3 16.61       4.58   0.02   0.96 -4.86
rain_days       3 16.72       4.68   0.02   0.97 -4.92
minT_maxW       4 17.03       5.00   0.01   0.99 -3.75
minT_maxR       4 18.24       6.20   0.01   1.00 -4.35
minT_x_meanWind 5 20.27       8.23   0.00   1.00 -3.93
minT_x_rainDays 5 21.95       9.92   0.00   1.00 -4.78
# AIC 2: with year
cand_effort_with_year <- list(
  year               = m_year,

  mean_temp          = m_mean_temp_year,
  max_temp           = m_max_temp_year,
  min_temp           = m_min_temp_year,
  total_rain         = m_total_rain_year,
  rain_days          = m_rain_days_year,
  max_rain           = m_max_rain_year,
  lag_total_rain     = m_lag_total_rain_year,
  mean_wind          = m_mean_wind_year,
  max_wind           = m_max_wind_year,
  hail_days          = m_hail_days_year,

  minT_totalR        = m_minT_totalR_year,
  minT_maxR          = m_minT_maxR_year,
  minT_maxW          = m_minT_maxW_year,
  hail_maxW          = m_hail_maxW_year,
  maxR_maxW          = m_maxR_maxW_year,
  maxT_rainDays      = m_maxT_rainDays_year,
  hail_rainDays      = m_hail_rainDays_year,

  minT_x_rainDays    = m_minT_x_rainDays_year,
  minT_x_meanWind    = m_minT_x_meanWind_year,
  lagR_x_maxT        = m_lagR_x_maxT_year
)

aic_effort_with_year <- aictab(cand.set = cand_effort_with_year, modnames = names(cand_effort_with_year))
aic_effort_with_year

Model selection based on AICc:

                K   AICc Delta_AICc AICcWt Cum.Wt    LL
max_rain        4 -29.07       0.00   0.19   0.19 19.31
year            3 -28.43       0.64   0.14   0.33 17.66
total_rain      4 -27.49       1.58   0.09   0.42 18.51
minT_maxR       5 -27.44       1.63   0.09   0.51 19.92
min_temp        4 -26.68       2.39   0.06   0.56 18.11
mean_temp       4 -26.59       2.48   0.06   0.62 18.06
max_temp        4 -26.42       2.65   0.05   0.67 17.98
maxR_maxW       5 -26.41       2.66   0.05   0.72 19.40
max_wind        4 -25.90       3.18   0.04   0.76 17.72
hail_days       4 -25.88       3.19   0.04   0.80 17.71
mean_wind       4 -25.87       3.21   0.04   0.84 17.70
lag_total_rain  4 -25.81       3.27   0.04   0.88 17.67
rain_days       4 -25.78       3.29   0.04   0.92 17.66
minT_totalR     5 -25.04       4.03   0.03   0.94 18.72
minT_maxW       5 -23.83       5.24   0.01   0.96 18.11
maxT_rainDays   5 -23.78       5.29   0.01   0.97 18.09
hail_maxW       5 -23.12       5.95   0.01   0.98 17.76
hail_rainDays   5 -23.03       6.04   0.01   0.99 17.72
lagR_x_maxT     6 -21.36       7.71   0.00   0.99 18.43
minT_x_meanWind 6 -21.30       7.77   0.00   1.00 18.40
minT_x_rainDays 6 -21.04       8.03   0.00   1.00 18.27

Check residuals

# identify top-ranked model (no year)
best_no_year_name <- aic_effort_no_year$Modnames[1]
best_no_year_E <- cand_effort_no_year[[best_no_year_name]]

par(mfrow = c(2,2))
plot(best_no_year_E)

# identify top-ranked model (with year)
best_with_year_name <- aic_effort_with_year$Modnames[1]
best_with_year_E <- cand_effort_with_year[[best_with_year_name]]

par(mfrow = c(2,2))
plot(best_with_year_E)

# year-only baseline and best with-year model
summary(m_year)$adj.r.squared
[1] 0.758928
summary(best_with_year_E)$adj.r.squared
[1] 0.7754688
# best no-year model
summary(best_no_year_E)$adj.r.squared
[1] 0.1105514

Figures

ggplot(effort_weather, aes(x = year, y = active_nests)) +
  geom_point(size = 2, colour = "grey30") +
  geom_smooth(method = "lm", se = TRUE, colour = "black") +
  scale_x_continuous(breaks = seq(min(effort_weather$year), max(effort_weather$year), by = 5)) +
  labs(
    x = "Year",
    y = "Number of active nests",
    title = "Annual breeding effort at Dronfield Nature Reserve"
  ) +
  theme_bw() +
  theme(panel.grid = element_blank())
`geom_smooth()` using formula = 'y ~ x'

prep_aic_plot <- function(aic_obj, set_label){
  df <- as.data.frame(aic_obj)

  # model names
  if ("Modnames" %in% names(df)) df$model <- df$Modnames
  if (!"model" %in% names(df)) df$model <- rownames(df)

  df |>
    dplyr::mutate(set = set_label) |>
    dplyr::rename(delta = Delta_AICc, weight = AICcWt) |>
    dplyr::select(set, model, delta, weight)
}

aic_plot_df <- dplyr::bind_rows(
  prep_aic_plot(aic_effort_with_year, "With year"),
  prep_aic_plot(aic_effort_no_year, "Without year")
) |>
  dplyr::filter(delta <= 3) |>
  dplyr::group_by(set) |>
  dplyr::mutate(model = forcats::fct_reorder(model, -delta)) |>
  dplyr::ungroup()

ggplot(aic_plot_df, aes(x = delta, y = model)) +
  geom_vline(xintercept = 2, linetype = 2) +
  geom_point(aes(size = weight), colour = "grey20") +
  facet_wrap(~set, scales = "free_y") +
  scale_size_continuous(range = c(2, 8)) +
  labs(
    x = expression(Delta*"AICc (lower is better)"),
    y = NULL,
    title = "Support for candidate models with and without year",
    subtitle = "Only models with ΔAICc ≤ 3 are shown; point size indicates AICc weight"
  ) +
  theme_bw() +
  theme(panel.grid = element_blank())

# Table 1(weather-only) 
aic_df <- as.data.frame(aic_effort_no_year)

# model ids
if ("Modnames" %in% names(aic_df)) {
  aic_df <- aic_df |> dplyr::mutate(model_id = Modnames)
} else {
  aic_df <- tibble::rownames_to_column(aic_df, "model_id")
}

table1 <- aic_df |>
  dplyr::rename(
    K = K,
    AICc = AICc,
    `ΔAICc` = Delta_AICc,
    w = AICcWt
  ) |>
  dplyr::select(model_id, K, AICc, `ΔAICc`, w) |>
  dplyr::arrange(`ΔAICc`) 

# add Adj R2 + slope(SE)
get_stats <- function(mod){
  adjr2 <- summary(mod)$adj.r.squared
  if (length(coef(mod)) == 2) {
    cf <- summary(mod)$coefficients
    slope_se <- paste0(round(cf[2,1], 3), " (", round(cf[2,2], 3), ")")
  } else {
    slope_se <- ""
  }
  tibble::tibble(`Adj. R²` = adjr2, `slope (SE)` = slope_se)
}

stats_df <- lapply(table1$model_id, function(mn) get_stats(cand_effort_no_year[[mn]])) |>
  dplyr::bind_rows() |>
  dplyr::mutate(model_id = table1$model_id)

table1 <- table1 |>
  dplyr::left_join(stats_df, by = "model_id") |>
  dplyr::mutate(
    AICc = round(AICc, 2),
    `ΔAICc` = round(`ΔAICc`, 2),
    w = round(w, 3),
    `Adj. R²` = round(`Adj. R²`, 3)
  )

# grouping adn neat names 
group_order <- c("Baseline","Rain","Temperature","Wind","Hail","Joint effects","Interactions")

model_group <- function(m){
  dplyr::case_when(
    m == "null" ~ "Baseline",
    m %in% c("total_rain","rain_days","max_rain","lag_total_rain") ~ "Rain",
    m %in% c("mean_temp","max_temp","min_temp") ~ "Temperature",
    m %in% c("mean_wind","max_wind") ~ "Wind",
    m %in% c("hail_days") ~ "Hail",
    m %in% c("minT_totalR","minT_maxR","minT_maxW","hail_maxW","maxR_maxW","maxT_rainDays","hail_rainDays") ~ "Joint effects",
    m %in% c("minT_x_rainDays","minT_x_meanWind","lagR_x_maxT") ~ "Interactions",
    TRUE ~ "Joint effects"
  )
}

pretty_model <- function(x){
  dplyr::recode(
    x,
    "null"           = "Intercept only",
    "total_rain"     = "Total rainfall",
    "rain_days"      = "Rainy days",
    "max_rain"       = "Maximum daily rainfall",
    "lag_total_rain" = "Previous year's total rainfall",
    "mean_temp"      = "Mean temperature",
    "max_temp"       = "Maximum temperature",
    "min_temp"       = "Minimum temperature",
    "mean_wind"      = "Mean wind speed",
    "max_wind"       = "Maximum wind speed",
    "hail_days"      = "Hail days",
    "minT_totalR"     = "Min temperature + total rainfall",
    "minT_maxR"       = "Min temperature + max daily rainfall",
    "minT_maxW"       = "Min temperature + max wind speed",
    "hail_maxW"       = "Hail days + max wind speed",
    "maxR_maxW"       = "Max daily rainfall + max wind speed",
    "maxT_rainDays"   = "Max temperature + rainy days",
    "hail_rainDays"   = "Hail days + rainy days",
    "minT_x_rainDays" = "Min temperature × rainy days",
    "minT_x_meanWind" = "Min temperature × mean wind speed",
    "lagR_x_maxT"     = "Prev. rainfall × max temperature",
    .default = x
  )
}

table1_clean <- table1 |>
  dplyr::mutate(
    Group = model_group(model_id),
    Model = pretty_model(model_id),
    Group = factor(Group, levels = group_order)
  ) |>
  dplyr::arrange(Group, `ΔAICc`) |>
  dplyr::select(Group, Model, K, AICc, `ΔAICc`, w, `Adj. R²`, `slope (SE)`)

#print 
group_sizes <- table1_clean |>
  dplyr::count(Group, name = "n") |>
  dplyr::mutate(
    start = cumsum(dplyr::lag(n, default = 0)) + 1,
    end   = cumsum(n)
  )

tbl <- knitr::kable(
  table1_clean |> dplyr::select(-Group),
  caption = "Table 1. Competitive weather-only candidate models (ΔAICc ≤ 4) explaining annual breeding effort (log number of active nests) at Dronfield. K = number of parameters; AICc = small-sample AIC; ΔAICc = difference from the best model; w = Akaike weight. Adjusted R² is shown for all models; slope (SE) is reported only for single-predictor models.",
  align = "l"
) |>
  kableExtra::kable_styling(full_width = FALSE)

for(i in seq_len(nrow(group_sizes))){
  tbl <- tbl |>
    kableExtra::pack_rows(as.character(group_sizes$Group[i]),
                          group_sizes$start[i],
                          group_sizes$end[i])
}

tbl
Table 1. Competitive weather-only candidate models (ΔAICc ≤ 4) explaining annual breeding effort (log number of active nests) at Dronfield. K = number of parameters; AICc = small-sample AIC; ΔAICc = difference from the best model; w = Akaike weight. Adjusted R² is shown for all models; slope (SE) is reported only for single-predictor models.
Model K AICc ΔAICc w Adj. R² slope (SE)
Baseline
Intercept only 2 14.26 2.22 0.057 0.000
Rain
Total rainfall 3 12.83 0.79 0.117 0.088 -0.001 (0)
Previous year's total rainfall 3 15.07 3.03 0.038 0.019 0 (0)
Maximum daily rainfall 3 15.74 3.70 0.027 -0.002 -0.003 (0.003)
Rainy days 3 16.72 4.68 0.017 -0.034 0 (0.003)
Temperature
Maximum temperature 3 12.77 0.74 0.120 0.089 0.046 (0.023)
Mean temperature 3 16.46 4.42 0.019 -0.026 0.031 (0.062)
Minimum temperature 3 16.61 4.58 0.018 -0.031 0.011 (0.036)
Wind
Maximum wind speed 3 14.39 2.35 0.054 0.041 -0.016 (0.011)
Mean wind speed 3 14.84 2.81 0.043 0.026 -0.18 (0.134)
Hail
Hail days 3 12.04 0.00 0.174 0.111 0.076 (0.035)
Joint effects
Hail days + max wind speed 4 13.16 1.12 0.099 0.123
Max temperature + rainy days 4 14.37 2.33 0.054 0.088
Hail days + rainy days 4 14.44 2.40 0.052 0.086
Min temperature + total rainfall 4 15.45 3.41 0.032 0.056
Max daily rainfall + max wind speed 4 15.84 3.81 0.026 0.044
Min temperature + max wind speed 4 17.03 5.00 0.014 0.006
Min temperature + max daily rainfall 4 18.24 6.20 0.008 -0.033
Interactions
Prev. rainfall × max temperature 5 15.71 3.67 0.028 0.100
Min temperature × mean wind speed 5 20.27 8.23 0.003 -0.043
Min temperature × rainy days 5 21.95 9.92 0.001 -0.101
digits_aic   <- 2
digits_r2    <- 3
digits_slope <- 3

group_order <- c("Baseline", "Rain", "Temperature", "Wind", "Hail", "Joint effects", "Interactions")

model_group_with_year <- function(m){
  dplyr::case_when(
    m %in% c("year") ~ "Baseline",
    m %in% c("total_rain","rain_days","max_rain","lag_total_rain") ~ "Rain",
    m %in% c("mean_temp","max_temp","min_temp") ~ "Temperature",
    m %in% c("mean_wind","max_wind") ~ "Wind",
    m %in% c("hail_days") ~ "Hail",
    m %in% c("minT_totalR","minT_maxR","minT_maxW","hail_maxW","maxR_maxW","maxT_rainDays","hail_rainDays") ~ "Joint effects",
    m %in% c("minT_x_rainDays","minT_x_meanWind","lagR_x_maxT") ~ "Interactions",
    TRUE ~ "Joint effects"
  )
}

pretty_model_with_year <- function(x){
  dplyr::recode(
    x,
    "year"           = "Year (baseline)",
    "total_rain"     = "Total rainfall",
    "rain_days"      = "Rainy days",
    "max_rain"       = "Maximum daily rainfall",
    "lag_total_rain" = "Previous year's total rainfall",
    "mean_temp"      = "Mean temperature",
    "max_temp"       = "Maximum temperature",
    "min_temp"       = "Minimum temperature",
    "mean_wind"      = "Mean wind speed",
    "max_wind"       = "Maximum wind speed",
    "hail_days"      = "Hail days",
    "minT_totalR"     = "Min temperature + total rainfall",
    "minT_maxR"       = "Min temperature + max daily rainfall",
    "minT_maxW"       = "Min temperature + max wind speed",
    "hail_maxW"       = "Hail days + max wind speed",
    "maxR_maxW"       = "Max daily rainfall + max wind speed",
    "maxT_rainDays"   = "Max temperature + rainy days",
    "hail_rainDays"   = "Hail days + rainy days",
    "minT_x_rainDays" = "Min temperature × rainy days",
    "minT_x_meanWind" = "Min temperature × mean wind speed",
    "lagR_x_maxT"     = "Prev. rainfall × max temperature",
    .default = x
  )
}

aic_df2 <- as.data.frame(aic_effort_with_year)

if ("Modnames" %in% names(aic_df2)) {
  aic_df2 <- aic_df2 |> dplyr::mutate(model = Modnames)
} else {
  aic_df2 <- tibble::rownames_to_column(aic_df2, "model")
}

aic_df2 <- aic_df2 |>
  dplyr::rename(
    K     = K,
    AICc  = AICc,
    dAICc = Delta_AICc,
    w     = AICcWt
  ) |>
  dplyr::select(model, K, AICc, dAICc, w) |>
  dplyr::arrange(dAICc)

get_stats_with_year <- function(mod){
  adjr2 <- summary(mod)$adj.r.squared

  if (length(coef(mod)) == 3) {   # intercept + year + ONE weather predictor
    cf <- summary(mod)$coefficients
    slope <- unname(cf[3, 1])
    se    <- unname(cf[3, 2])
    slope_se <- paste0(
      round(slope, digits_slope),
      " (", round(se, digits_slope), ")"
    )
  } else {
    slope_se <- ""
  }

  tibble::tibble(adj_R2 = adjr2, slope_SE = slope_se)
}


stats_df2 <- lapply(aic_df2$model, function(mn) get_stats_with_year(cand_effort_with_year[[mn]])) |>
  dplyr::bind_rows() |>
  dplyr::mutate(model = aic_df2$model)

table2 <- aic_df2 |>
  dplyr::left_join(stats_df2, by = "model") |>
  dplyr::mutate(
    AICc   = round(AICc, digits_aic),
    dAICc  = round(dAICc, digits_aic),
    w      = round(w, 3),
    adj_R2 = round(adj_R2, digits_r2)
  ) |>
  dplyr::mutate(
    model_id = model,
    Group = model_group_with_year(model_id),
    Model = pretty_model_with_year(model_id),
    Group = factor(Group, levels = group_order)
  ) |>
  dplyr::arrange(Group, dAICc) |>
  dplyr::rename(
    `ΔAICc` = dAICc,
    `Adj. R²` = adj_R2,
    `slope (SE)` = slope_SE
  ) |>
  dplyr::select(Group, Model, K, AICc, `ΔAICc`, w, `Adj. R²`, `slope (SE)`)

group_sizes2 <- table2 |>
  dplyr::count(Group, name = "n") |>
  dplyr::filter(n > 0) |>
  dplyr::mutate(
    start = cumsum(dplyr::lag(n, default = 0)) + 1,
    end   = cumsum(n)
  )

tbl2 <- knitr::kable(
  table2 |> dplyr::select(-Group),
  caption = "Table 2. Candidate models including year explaining annual breeding effort (log number of active nests) at Dronfield. K = number of parameters; AICc = small-sample AIC; ΔAICc = difference from the best model; w = Akaike weight. Adjusted R² is shown for all models; slope (SE) is reported for each model (non-intercept terms shown as estimate (SE)).",
  align = "l"
) |>
  kableExtra::kable_styling(full_width = FALSE)

for(i in seq_len(nrow(group_sizes2))){
  tbl2 <- tbl2 |>
    kableExtra::pack_rows(
      as.character(group_sizes2$Group[i]),
      group_sizes2$start[i],
      group_sizes2$end[i]
    )
}

tbl2
Table 2. Candidate models including year explaining annual breeding effort (log number of active nests) at Dronfield. K = number of parameters; AICc = small-sample AIC; ΔAICc = difference from the best model; w = Akaike weight. Adjusted R² is shown for all models; slope (SE) is reported for each model (non-intercept terms shown as estimate (SE)).
Model K AICc ΔAICc w Adj. R² slope (SE)
Baseline
Year (baseline) 3 -28.43 0.64 0.140 0.759
Rain
Maximum daily rainfall 4 -29.07 0.00 0.193 0.775 -0.002 (0.001)
Total rainfall 4 -27.49 1.58 0.088 0.764 0 (0)
Previous year's total rainfall 4 -25.81 3.27 0.038 0.751 0 (0)
Rainy days 4 -25.78 3.29 0.037 0.750 0 (0.002)
Temperature
Minimum temperature 4 -26.68 2.39 0.058 0.757 0.016 (0.017)
Mean temperature 4 -26.59 2.48 0.056 0.757 0.026 (0.03)
Maximum temperature 4 -26.42 2.65 0.051 0.755 0.01 (0.013)
Wind
Maximum wind speed 4 -25.90 3.18 0.039 0.751 -0.002 (0.006)
Mean wind speed 4 -25.87 3.21 0.039 0.751 -0.019 (0.07)
Hail
Hail days 4 -25.88 3.19 0.039 0.751 0.006 (0.02)
Joint effects
Min temperature + max daily rainfall 5 -27.44 1.63 0.085 0.776
Max daily rainfall + max wind speed 5 -26.41 2.66 0.051 0.769
Min temperature + total rainfall 5 -25.04 4.03 0.026 0.758
Min temperature + max wind speed 5 -23.83 5.24 0.014 0.749
Max temperature + rainy days 5 -23.78 5.29 0.014 0.748
Hail days + max wind speed 5 -23.12 5.95 0.010 0.743
Hail days + rainy days 5 -23.03 6.04 0.009 0.742
Interactions
Prev. rainfall × max temperature 6 -21.36 7.71 0.004 0.744
Min temperature × mean wind speed 6 -21.30 7.77 0.004 0.744
Min temperature × rainy days 6 -21.04 8.03 0.003 0.741

Objective 2: Breeding success

objective2_table <- tribble(
  ~`Hypothesis / justification`, ~Model,

  # ---- Temperature ----
  "**Temp**", "",
  "Extreme heat events cause chick heat stress and dehydration, reducing fledging success",
  "success ~ max_temp + (1 | year)",

  "Cold extremes increase thermoregulatory demand in chicks, increasing mortality risk",
  "success ~ min_temp + (1 | year)",

  # ---- Rain ----
  "**Rain**", "",
  "Prolonged rainfall reduces adult foraging efficiency, limiting food delivery to chicks and reducing breeding success",
  "success ~ rain_days + (1 | year)",

  "Intense rainfall events cause acute nest disturbance and chick exposure, reducing breeding success",
  "success ~ max_rain + (1 | year)",

  "Overall wet years reduce provisioning efficiency across the breeding season, lowering breeding success",
  "success ~ total_rain + (1 | year)",

  # ---- Wind ----
  "**Wind**", "",
  "Extreme winds increase nest exposure and disrupt provisioning, reducing chick survival",
  "success ~ max_wind + (1 | year)",

  # ---- Hail ----
  "**Hail**", "",
  "Severe hail events increase nest destruction and chick mortality, reducing breeding success",
  "success ~ hail_days + (1 | year)",

  # ---- Joint Effects ----
  "**Joint effects**", "",
  "Severe storm exposure increases nest destruction and chick mortality, reducing breeding success",
  "success ~ hail_days + max_wind + (1 | year)",

  "Prolonged wet conditions combined with storm events increase chick exposure and limit provisioning, reducing breeding success",
  "success ~ rain_days + hail_days + (1 | year)",

  "Heat stress combined with prolonged wet conditions reduces chick thermoregulation and food delivery, lowering breeding success",
  "success ~ max_temp + rain_days + (1 | year)",

  # ---- Interactive Effects ----
  "**Interactive effects**", "",
  "The negative effect of heat stress on breeding success is amplified during persistently wet conditions that limit adult foraging",
  "success ~ max_temp * rain_days + (1 | year)",

  "Cold stress on chicks is amplified during extreme wind events due to increased exposure and heat loss, reducing breeding success",
  "success ~ min_temp * max_wind + (1 | year)"
)

kable(
  objective2_table,
  align = c("l", "l"),
  caption = "Objective 2 candidate model set for breeding success. All models are fitted as binomial GLMMs: cbind(success_nests, failed_nests) with a random intercept for year (1 | year)"
)
Objective 2 candidate model set for breeding success. All models are fitted as binomial GLMMs: cbind(success_nests, failed_nests) with a random intercept for year (1 | year)
Hypothesis / justification Model
Temp
Extreme heat events cause chick heat stress and dehydration, reducing fledging success success ~ max_temp + (1 | year)
Cold extremes increase thermoregulatory demand in chicks, increasing mortality risk success ~ min_temp + (1 | year)
Rain
Prolonged rainfall reduces adult foraging efficiency, limiting food delivery to chicks and reducing breeding success success ~ rain_days + (1 | year)
Intense rainfall events cause acute nest disturbance and chick exposure, reducing breeding success success ~ max_rain + (1 | year)
Overall wet years reduce provisioning efficiency across the breeding season, lowering breeding success success ~ total_rain + (1 | year)
Wind
Extreme winds increase nest exposure and disrupt provisioning, reducing chick survival success ~ max_wind + (1 | year)
Hail
Severe hail events increase nest destruction and chick mortality, reducing breeding success success ~ hail_days + (1 | year)
Joint effects
Severe storm exposure increases nest destruction and chick mortality, reducing breeding success success ~ hail_days + max_wind + (1 | year)
Prolonged wet conditions combined with storm events increase chick exposure and limit provisioning, reducing breeding success success ~ rain_days + hail_days + (1 | year)
Heat stress combined with prolonged wet conditions reduces chick thermoregulation and food delivery, lowering breeding success success ~ max_temp + rain_days + (1 | year)
Interactive effects
The negative effect of heat stress on breeding success is amplified during persistently wet conditions that limit adult foraging success ~ max_temp * rain_days + (1 | year)
Cold stress on chicks is amplified during extreme wind events due to increased exposure and heat loss, reducing breeding success success ~ min_temp * max_wind + (1 | year)

success ~ 1 + (1 | year)

library(lme4)

success_df <- success_df |> mutate(year = factor(year))

m0 <- glmer(
  cbind(success_nests, failed_nests) ~ 1 + (1 | year),
  family = binomial,
  data = success_df
)
summary(m0)
Generalized linear mixed model fit by maximum likelihood (Laplace
  Approximation) [glmerMod]
 Family: binomial  ( logit )
Formula: cbind(success_nests, failed_nests) ~ 1 + (1 | year)
   Data: success_df

      AIC       BIC    logLik -2*log(L)  df.resid 
    229.0     232.0    -112.5     225.0        30 

Scaled residuals: 
     Min       1Q   Median       3Q      Max 
-1.25016 -0.28770  0.00005  0.48287  1.24098 

Random effects:
 Groups Name        Variance Std.Dev.
 year   (Intercept) 0.1292   0.3594  
Number of obs: 32, groups:  year, 32

Fixed effects:
            Estimate Std. Error z value Pr(>|z|)   
(Intercept)  0.21944    0.07595   2.889  0.00386 **
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
library(lme4)

success_df <- success_df |> dplyr::mutate(year = factor(year))

# Baseline (null) model

m0_success <- glmer(
  cbind(success_nests, failed_nests) ~ 1 + (1 | year),
  family = binomial,
  data = success_df
)

# Single predictor models

m_succ_maxT <- glmer(
  cbind(success_nests, failed_nests) ~ z_max_temp + (1 | year),
  family = binomial, data = success_df
)

m_succ_minT <- glmer(
  cbind(success_nests, failed_nests) ~ z_min_temp + (1 | year),
  family = binomial, data = success_df
)

m_succ_rainD <- glmer(
  cbind(success_nests, failed_nests) ~ z_rain_days + (1 | year),
  family = binomial, data = success_df
)

m_succ_maxR <- glmer(
  cbind(success_nests, failed_nests) ~ z_max_rain + (1 | year),
  family = binomial, data = success_df
)

m_succ_totR <- glmer(
  cbind(success_nests, failed_nests) ~ z_total_rain + (1 | year),
  family = binomial, data = success_df
)

m_succ_maxW <- glmer(
  cbind(success_nests, failed_nests) ~ z_max_wind + (1 | year),
  family = binomial, data = success_df
)

m_succ_hail <- glmer( 
  cbind(success_nests, failed_nests) ~ z_hail_days + (1 | year),
  family = binomial, data = success_df
)

# Joint effects

m_succ_hailW <- glmer(
  cbind(success_nests, failed_nests) ~ z_hail_days + z_max_wind + (1 | year),
  family = binomial, data = success_df
)

m_succ_rainH <- glmer(
  cbind(success_nests, failed_nests) ~ z_rain_days + z_hail_days + (1 | year),
  family = binomial, data = success_df
)

m_succ_heatR <- glmer(
  cbind(success_nests, failed_nests) ~ z_max_temp + z_rain_days + (1 | year),
  family = binomial, data = success_df
)

# Interactions

m_succ_int_heat_rain <- glmer(
  cbind(success_nests, failed_nests) ~ z_max_temp * z_rain_days + (1 | year),
  family = binomial, data = success_df
)

m_succ_int_cold_wind <- glmer(
  cbind(success_nests, failed_nests) ~ z_min_temp * z_max_wind + (1 | year),
  family = binomial, data = success_df
)
# Candidate set + AIC

library(AICcmodavg)

cand_success <- list(
  null_re      = m0_success,

  max_temp     = m_succ_maxT,
  min_temp     = m_succ_minT,
  rain_days    = m_succ_rainD,
  max_rain     = m_succ_maxR,
  total_rain   = m_succ_totR,
  max_wind     = m_succ_maxW,
  hail_days    = m_succ_hail,

  hail_maxW    = m_succ_hailW,
  rain_hail    = m_succ_rainH,
  maxT_rainD   = m_succ_heatR,

  int_heatRain = m_succ_int_heat_rain,
  int_coldWind = m_succ_int_cold_wind
)

aic_success <- aictab(
  cand.set  = cand_success,
  modnames  = names(cand_success)
)

aic_success

Model selection based on AICc:

             K   AICc Delta_AICc AICcWt Cum.Wt      LL
min_temp     3 226.93       0.00   0.43   0.43 -110.03
null_re      2 229.45       2.52   0.12   0.56 -112.52
hail_days    3 230.45       3.53   0.07   0.63 -111.80
max_temp     3 230.85       3.92   0.06   0.69 -111.99
max_rain     3 231.01       4.09   0.06   0.75 -112.08
rain_days    3 231.61       4.68   0.04   0.79 -112.38
total_rain   3 231.65       4.72   0.04   0.83 -112.40
max_wind     3 231.73       4.81   0.04   0.87 -112.44
int_coldWind 5 232.04       5.11   0.03   0.91 -109.87
maxT_rainD   4 232.10       5.18   0.03   0.94 -111.31
rain_hail    4 232.30       5.38   0.03   0.97 -111.41
hail_maxW    4 232.67       5.74   0.02   0.99 -111.59
int_heatRain 5 234.93       8.01   0.01   1.00 -111.31
vc_year <- function(mod) {
  vc <- as.data.frame(VarCorr(mod))
  vc$vcov[vc$grp == "year"][1]
}

v0     <- vc_year(m0_success)
v_minT <- vc_year(m_succ_minT)

prop_between_year_explained <- (v0 - v_minT) / v0
prop_between_year_explained
[1] 0.2098479

Figures and tables

library(dplyr)
library(tibble)
library(purrr)
library(knitr)
library(kableExtra)

Attaching package: 'kableExtra'
The following object is masked from 'package:dplyr':

    group_rows
library(lme4)

aic_df <- as.data.frame(aic_success)

if ("Modnames" %in% names(aic_df)) {
  aic_df <- aic_df |> mutate(model_id = Modnames)
} else {
  aic_df <- tibble::rownames_to_column(aic_df, "model_id")
}

table_success <- aic_df |>
  rename(
    K = K,
    AICc = AICc,
    `ΔAICc` = Delta_AICc,
    w = AICcWt
  ) |>
  dplyr::select(model_id, K, AICc, `ΔAICc`, w) |>
  arrange(`ΔAICc`)

vc_year <- function(mod) {
  vc <- as.data.frame(VarCorr(mod))
  vc$vcov[vc$grp == "year"][1]
}

v0 <- vc_year(m0_success)

var_expl_tbl <- tibble(
  model_id = names(cand_success),
  `Variance explained` = map_dbl(cand_success, ~ {
    v <- vc_year(.x)
    (v0 - v) / v0
  })
) |>
  mutate(`Variance explained` = round(`Variance explained`, 3))

table_success <- table_success |>
  left_join(var_expl_tbl, by = "model_id")

get_slope_se <- function(mod) {
  cf <- summary(mod)$coefficients
  if (nrow(cf) == 2) {
    paste0(round(cf[2, 1], 3), " (", round(cf[2, 2], 3), ")")
  } else {
    ""
  }
}

slope_tbl <- tibble(
  model_id = table_success$model_id,
  `slope (SE)` = map_chr(table_success$model_id, ~ get_slope_se(cand_success[[.x]]))
)

table_success <- table_success |>
  left_join(slope_tbl, by = "model_id") |>
  mutate(
    AICc = round(AICc, 2),
    `ΔAICc` = round(`ΔAICc`, 2),
    w = round(w, 3)
  )

group_order <- c("Baseline", "Temperature", "Rain", "Wind", "Hail", "Joint effects", "Interactions")

model_group_succ <- function(m) {
  case_when(
    m %in% c("null_re", "null", "null_success") ~ "Baseline",
    m %in% c("min_temp", "max_temp") ~ "Temperature",
    m %in% c("rain_days", "max_rain", "total_rain") ~ "Rain",
    m %in% c("max_wind") ~ "Wind",
    m %in% c("hail_days") ~ "Hail",
    m %in% c("hail_maxW", "rain_hail", "maxT_rainD") ~ "Joint effects",
    m %in% c("int_heatRain", "int_coldWind") ~ "Interactions",
    TRUE ~ "Joint effects"
  )
}

neat_succ <- function(x) {
  recode(
    x,
    "null_re"      = "Intercept only",
    "min_temp"     = "Minimum temperature",
    "max_temp"     = "Maximum temperature",
    "rain_days"    = "Rainy days",
    "max_rain"     = "Maximum daily rainfall",
    "total_rain"   = "Total rainfall",
    "max_wind"     = "Maximum wind speed",
    "hail_days"    = "Hail days",
    "hail_maxW"    = "Hail days + max wind speed",
    "rain_hail"    = "Rainy days + hail days",
    "maxT_rainD"   = "Max temperature + rainy days",
    "int_heatRain" = "Max temperature × rainy days",
    "int_coldWind" = "Min temperature × max wind speed",
    .default = x
  )
}

table_success_clean <- table_success |>
  mutate(
    Group = model_group_succ(model_id),
    Model = neat_succ(model_id),
    Group = factor(Group, levels = group_order)
  ) |>
  arrange(Group, `ΔAICc`) |>
  dplyr::select(Group, Model, K, AICc, `ΔAICc`, w, `Variance explained`, `slope (SE)`)

group_sizes <- table_success_clean |>
  count(Group, name = "n") |>
  filter(n > 0) |>
  mutate(
    start = cumsum(lag(n, default = 0)) + 1,
    end = cumsum(n)
  )

tbl_succ <- kable(
  table_success_clean |> dplyr::select(-Group),
  caption = "Model selection results for breeding success (binomial GLMMs with a random intercept for year). Variance explained is the proportional reduction in the estimated between-year variance relative to the null model. Slope (SE) is shown only for single-predictor models.",
  align = "l"
) |>
  kable_styling(full_width = FALSE)

for (i in seq_len(nrow(group_sizes))) {
  tbl_succ <- tbl_succ |>
    pack_rows(
      as.character(group_sizes$Group[i]),
      group_sizes$start[i],
      group_sizes$end[i]
    )
}

tbl_succ
Model selection results for breeding success (binomial GLMMs with a random intercept for year). Variance explained is the proportional reduction in the estimated between-year variance relative to the null model. Slope (SE) is shown only for single-predictor models.
Model K AICc ΔAICc w Variance explained slope (SE)
Baseline
Intercept only 2 229.45 2.52 0.123 0.000
Temperature
Minimum temperature 3 226.93 0.00 0.435 0.210 0.167 (0.072)
Maximum temperature 3 230.85 3.92 0.061 0.042 -0.077 (0.075)
Rain
Maximum daily rainfall 3 231.01 4.09 0.056 0.034 0.075 (0.079)
Rainy days 3 231.61 4.68 0.042 0.014 -0.041 (0.076)
Total rainfall 3 231.65 4.72 0.041 0.008 0.038 (0.076)
Wind
Maximum wind speed 3 231.73 4.81 0.039 0.003 -0.031 (0.077)
Hail
Hail days 3 230.45 3.53 0.074 0.052 -0.088 (0.073)
Joint effects
Max temperature + rainy days 4 232.10 5.18 0.033 0.104
Rainy days + hail days 4 232.30 5.38 0.030 0.087
Hail days + max wind speed 4 232.67 5.74 0.025 0.062
Interactions
Min temperature × max wind speed 5 232.04 5.11 0.034 0.225
Max temperature × rainy days 5 234.93 8.01 0.008 0.104
minT_model <- m_succ_minT


success_df_plot <- success_df |>
  mutate(
    success_rate = success_nests / (success_nests + failed_nests),
    year_num = as.numeric(as.character(year))
  )

yrs <- sort(unique(success_df_plot$year_num[is.finite(success_df_plot$year_num)]))
breaks_5y <- yrs[seq(1, length(yrs), by = 5)]  #

ggplot(success_df_plot, aes(x = year_num, y = success_rate)) +
  geom_point(size = 2, colour = "grey30") +
  geom_smooth(
    method = "lm",
    se = TRUE,
    colour = "grey30",
    fill = "grey85"
  ) +
  scale_y_continuous(limits = c(0, 1)) +
  scale_x_continuous(breaks = breaks_5y) +
  labs(
    x = "Year",
    y = "Breeding success",
    title = "Annual breeding success (proportion of successful nests) of White-backed Vultures"
  ) +
  theme_bw() +
  theme(panel.grid = element_blank())
`geom_smooth()` using formula = 'y ~ x'

obs_df <- success_df |>
  mutate(
    success_rate = success_nests / (success_nests + failed_nests),
    se = sqrt(success_rate * (1 - success_rate) /
                (success_nests + failed_nests))
  )
minT_model <- m_succ_minT

pred_df <- success_df |>
  mutate(
    pred = predict(
      minT_model,
      newdata = success_df,
      type = "response",
      re.form = NA
    )
  )
obs_df <- success_df |>
  mutate(
    year_num = as.numeric(gsub("[^0-9]", "", as.character(year))),
    success_rate = success_nests / (success_nests + failed_nests),
    se = sqrt(success_rate * (1 - success_rate) / (success_nests + failed_nests))
  )

minT_model <- m_succ_minT

pred_df <- success_df |>
  mutate(
    year_num = as.numeric(gsub("[^0-9]", "", as.character(year))),
    pred = predict(minT_model, newdata = success_df, type = "response", re.form = NA)
  )


ggplot() +
  geom_point(
    data = obs_df,
    aes(x = year_num, y = success_rate),
    size = 2, shape = 21, fill = "white", colour = "grey30"
  ) +
  geom_errorbar(
    data = obs_df,
    aes(x = year_num,
        ymin = success_rate - se,
        ymax = success_rate + se),
    width = 0, colour = "grey30"
  ) +
  geom_line(
    data = pred_df,
    aes(x = year_num, y = pred),
    linewidth = 1, colour = "black"
  ) +
  scale_y_continuous(limits = c(0, 1)) +
  scale_x_continuous(breaks = seq(1995, 2025, by = 5)) +
  labs(
    x = "Year",
    y = "Breeding success",
    title = "Observed and temperature-predicted breeding success",
    subtitle = "Points show annual observed success (± SE); line shows predictions from the minimum temperature model"
  ) +
  theme_bw() +
  theme(panel.grid = element_blank())

figure shows that years with warmer minimum temperatures tend to have slightly higher breeding success, but temperature alone does not explain most of the year-to-year variation. Given the temperature in each year, this is what the model predicts breeding success should be


Objective 3

hist(timing_weather_60$lay_DOY)

qqnorm(timing_weather_60$lay_DOY)
qqline(timing_weather_60$lay_DOY)

objective3_table <- tribble(
  ~`Hypothesis / justification`, ~Model,

  # ---- Temperature (pre-laying cues) ----
  "**Temp (pre-laying)**", "",
  "Short-term warm conditions prior to laying may act as a cue for the initiation of breeding.",
  "lay_DOY ~ mean_temp_60",

  "Cold pre-laying conditions delay breeding due to increased energetic costs",
  "lay_DOY ~ min_temp_60",


  # ---- Rain (pre-laying cues) ----
  "**Rain (pre-laying)**", "",
  "Prolonged wet conditions prior to laying delay breeding by limiting foraging efficiency",
  "lay_DOY ~ rain_days_60",

  "Higher rainfall prior to laying signals improved food availability and advances breeding timing",
  "lay_DOY ~ total_rain_60",

  # ---- Wind (pre-laying flight conditions) ----
  "**Wind (pre-laying)**", "",
  "Persistent windy conditions prior to laying increase flight costs and delay breeding timing",
  "lay_DOY ~ mean_wind_60",

  # ---- Lagged effects (carry-over cues) ----
  "**Lagged effects (previous year)**", "",
  "Higher rainfall in the previous year improves food availability and adult condition entering the breeding season, allowing vultures to initiate breeding earlier",
  "lay_DOY ~ lag_total_rain",

  "Prolonged wet conditions in the previous year may influence adult condition and carry over to affect laying timing",
  "lay_DOY ~ lag_rain_days",

  "Thermal conditions in the previous year influence physiological state entering the breeding season",
  "lay_DOY ~ lag_mean_temp",

  # ---- Joint effects ----
  "**Joint effects**", "",
  "Temperature and rainfall cues prior to laying jointly influence breeding timing",
  "lay_DOY ~ mean_temp_60 + rain_days_60",

  "Carry-over effects from the previous year interact with current pre-laying rainfall cues",
  "lay_DOY ~ lag_total_rain + rain_days_60",

  "Previous-year rainfall and current temperature jointly influence breeding timing",
  "lay_DOY ~ lag_total_rain + mean_temp_60",

  # ---- Interactive effects ----
  "**Interactive effects**", "",
  "The effect of temperature cues on breeding timing depends on pre-laying rainfall conditions",
  "lay_DOY ~ mean_temp_60 * rain_days_60",

  "Carry-over effects from the previous year modify responses to current pre-laying temperature cues",
  "lay_DOY ~ lag_total_rain * mean_temp_60"
)

kable(
  objective3_table,
  align = c("l", "l"),
  caption = "Objective 3 candidate model set for breeding timing (laying date, day of year)."
)
Objective 3 candidate model set for breeding timing (laying date, day of year).
Hypothesis / justification Model
Temp (pre-laying)
Short-term warm conditions prior to laying may act as a cue for the initiation of breeding. lay_DOY ~ mean_temp_60
Cold pre-laying conditions delay breeding due to increased energetic costs lay_DOY ~ min_temp_60
Rain (pre-laying)
Prolonged wet conditions prior to laying delay breeding by limiting foraging efficiency lay_DOY ~ rain_days_60
Higher rainfall prior to laying signals improved food availability and advances breeding timing lay_DOY ~ total_rain_60
Wind (pre-laying)
Persistent windy conditions prior to laying increase flight costs and delay breeding timing lay_DOY ~ mean_wind_60
Lagged effects (previous year)
Higher rainfall in the previous year improves food availability and adult condition entering the breeding season, allowing vultures to initiate breeding earlier lay_DOY ~ lag_total_rain
Prolonged wet conditions in the previous year may influence adult condition and carry over to affect laying timing lay_DOY ~ lag_rain_days
Thermal conditions in the previous year influence physiological state entering the breeding season lay_DOY ~ lag_mean_temp
Joint effects
Temperature and rainfall cues prior to laying jointly influence breeding timing lay_DOY ~ mean_temp_60 + rain_days_60
Carry-over effects from the previous year interact with current pre-laying rainfall cues lay_DOY ~ lag_total_rain + rain_days_60
Previous-year rainfall and current temperature jointly influence breeding timing lay_DOY ~ lag_total_rain + mean_temp_60
Interactive effects
The effect of temperature cues on breeding timing depends on pre-laying rainfall conditions lay_DOY ~ mean_temp_60 * rain_days_60
Carry-over effects from the previous year modify responses to current pre-laying temperature cues lay_DOY ~ lag_total_rain * mean_temp_60
# 60 days prior
timing_weather_60 <- timing_weather_60 |>
  left_join(
    weather_yearly_lag |>
      dplyr::select(
        YEAR,
        lag_total_rain,
        lag_rain_days,
        lag_mean_temp
      ),
    by = c("year" = "YEAR")
  )

# 30 days prior
timing_weather_30 <- timing_weather_30 |>
  left_join(
    weather_yearly_lag |>
      dplyr::select(
        YEAR,
        lag_total_rain,
        lag_rain_days,
        lag_mean_temp
      ),
    by = c("year" = "YEAR")
  )
cand_timing_30 <- list(

  meanT_30   = lm(lay_DOY ~ mean_temp_30, data = timing_weather_30),
  maxT_30    = lm(lay_DOY ~ max_temp_30,  data = timing_weather_30),
  minT_30    = lm(lay_DOY ~ min_temp_30,  data = timing_weather_30),

  rainDays_30 = lm(lay_DOY ~ rain_days_30,  data = timing_weather_30),
  totalR_30   = lm(lay_DOY ~ total_rain_30, data = timing_weather_30),

  meanW_30 = lm(lay_DOY ~ mean_wind_30, data = timing_weather_30),
  maxW_30  = lm(lay_DOY ~ max_wind_30,  data = timing_weather_30),

  # lagged (previous year)
  lagTotR  = lm(lay_DOY ~ lag_total_rain, data = timing_weather_30),
  lagRainD = lm(lay_DOY ~ lag_rain_days,  data = timing_weather_30),
  lagMeanT = lm(lay_DOY ~ lag_mean_temp,  data = timing_weather_30),

  # additive
  temp_rain = lm(lay_DOY ~ mean_temp_30 + rain_days_30, data = timing_weather_30),
  temp_wind = lm(lay_DOY ~ mean_temp_30 + mean_wind_30, data = timing_weather_30),
  rain_wind = lm(lay_DOY ~ rain_days_30 + mean_wind_30, data = timing_weather_30),

  # interactions
  temp_x_rain = lm(lay_DOY ~ mean_temp_30 * rain_days_30, data = timing_weather_30),
  minT_x_wind = lm(lay_DOY ~ min_temp_30 * mean_wind_30, data = timing_weather_30),
  lagR_x_temp = lm(lay_DOY ~ lag_total_rain * mean_temp_30, data = timing_weather_30)
)
aic_timing_30 <- aictab(
  cand.set = cand_timing_30,
  modnames = names(cand_timing_30)
)
cand_timing_60 <- list(

  meanT_60   = lm(lay_DOY ~ mean_temp_60, data = timing_weather_60),
  maxT_60    = lm(lay_DOY ~ max_temp_60,  data = timing_weather_60),
  minT_60    = lm(lay_DOY ~ min_temp_60,  data = timing_weather_60),

  rainDays_60 = lm(lay_DOY ~ rain_days_60,  data = timing_weather_60),
  totalR_60   = lm(lay_DOY ~ total_rain_60, data = timing_weather_60),

  meanW_60 = lm(lay_DOY ~ mean_wind_60, data = timing_weather_60),
  maxW_60  = lm(lay_DOY ~ max_wind_60,  data = timing_weather_60),

  # lagged (previous year)
  lagTotR  = lm(lay_DOY ~ lag_total_rain, data = timing_weather_60),
  lagRainD = lm(lay_DOY ~ lag_rain_days,  data = timing_weather_60),
  lagMeanT = lm(lay_DOY ~ lag_mean_temp,  data = timing_weather_60),

  # additive
  temp_rain = lm(lay_DOY ~ mean_temp_60 + rain_days_60, data = timing_weather_60),
  temp_wind = lm(lay_DOY ~ mean_temp_60 + mean_wind_60, data = timing_weather_60),
  rain_wind = lm(lay_DOY ~ rain_days_60 + mean_wind_60, data = timing_weather_60),

  # interactions
  temp_x_rain = lm(lay_DOY ~ mean_temp_60 * rain_days_60, data = timing_weather_60),
  minT_x_wind = lm(lay_DOY ~ min_temp_60 * mean_wind_60, data = timing_weather_60),
  lagR_x_temp = lm(lay_DOY ~ lag_total_rain * mean_temp_60, data = timing_weather_60)
)
aic_timing_60 <- aictab(
  cand.set = cand_timing_60,
  modnames = names(cand_timing_60)
)
aic_timing_30

Model selection based on AICc:

            K     AICc Delta_AICc AICcWt Cum.Wt       LL
temp_x_rain 5 11336.00       0.00      1      1 -5662.98
temp_rain   4 11349.65      13.65      0      1 -5670.81
lagR_x_temp 5 11373.62      37.62      0      1 -5681.79
minT_x_wind 5 11438.51     102.51      0      1 -5714.24
minT_30     3 11529.94     193.94      0      1 -5761.96
temp_wind   4 11571.91     235.91      0      1 -5781.94
meanT_30    3 11600.33     264.33      0      1 -5797.16
maxT_30     3 11820.41     484.41      0      1 -5907.20
rain_wind   4 12224.48     888.48      0      1 -6108.23
lagTotR     3 12278.29     942.29      0      1 -6136.14
lagRainD    3 12289.62     953.61      0      1 -6141.80
lagMeanT    3 12292.70     956.70      0      1 -6143.34
rainDays_30 3 12360.96    1024.96      0      1 -6177.47
totalR_30   3 12400.54    1064.54      0      1 -6197.26
maxW_30     3 12506.89    1170.89      0      1 -6250.44
meanW_30    3 12516.32    1180.31      0      1 -6255.15
summary(cand_timing_30$temp_x_rain)$r.squared
[1] 0.5628433
aic_timing_60

Model selection based on AICc:

            K     AICc Delta_AICc AICcWt Cum.Wt       LL
temp_x_rain 5 10760.71       0.00    0.9    0.9 -5375.33
temp_rain   4 10765.17       4.46    0.1    1.0 -5378.57
lagR_x_temp 5 10889.01     128.31    0.0    1.0 -5439.49
minT_x_wind 5 10945.22     184.51    0.0    1.0 -5467.59
minT_60     3 10990.38     229.67    0.0    1.0 -5492.18
temp_wind   4 11100.18     339.48    0.0    1.0 -5546.08
meanT_60    3 11112.60     351.89    0.0    1.0 -5553.29
maxT_60     3 11394.57     633.86    0.0    1.0 -5694.28
rain_wind   4 12085.12    1324.41    0.0    1.0 -6038.55
rainDays_60 3 12211.31    1450.61    0.0    1.0 -6102.65
lagTotR     3 12278.29    1517.58    0.0    1.0 -6136.14
totalR_60   3 12282.33    1521.62    0.0    1.0 -6138.15
lagRainD    3 12289.62    1528.91    0.0    1.0 -6141.80
lagMeanT    3 12292.70    1531.99    0.0    1.0 -6143.34
maxW_60     3 12533.72    1773.01    0.0    1.0 -6263.85
meanW_60    3 12538.92    1778.21    0.0    1.0 -6266.45
summary(cand_timing_60$temp_x_rain)$r.squared
[1] 0.7052922
m_year_trend <- lm(lay_DOY ~ year, data = timing_weather_30)
summary(m_year_trend)

Call:
lm(formula = lay_DOY ~ year, data = timing_weather_30)

Residuals:
     Min       1Q   Median       3Q      Max 
-145.506  -10.576   -2.616    7.142  142.283 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)
(Intercept) 51.56562   99.65756   0.517    0.605
year         0.05049    0.04958   1.018    0.309

Residual standard error: 18.18 on 1458 degrees of freedom
Multiple R-squared:  0.0007109, Adjusted R-squared:  2.55e-05 
F-statistic: 1.037 on 1 and 1458 DF,  p-value: 0.3086
summary(m_year_trend)$r.squared
[1] 0.0007108822
m_year_trend60 <- lm(lay_DOY ~ year, data = timing_weather_60)
summary(m_year_trend60)

Call:
lm(formula = lay_DOY ~ year, data = timing_weather_60)

Residuals:
     Min       1Q   Median       3Q      Max 
-145.506  -10.576   -2.616    7.142  142.283 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)
(Intercept) 51.56562   99.65756   0.517    0.605
year         0.05049    0.04958   1.018    0.309

Residual standard error: 18.18 on 1458 degrees of freedom
Multiple R-squared:  0.0007109, Adjusted R-squared:  2.55e-05 
F-statistic: 1.037 on 1 and 1458 DF,  p-value: 0.3086

Figures and Tables

cand_timing_30 <- cand_timing_30[names(cand_timing_30) != "null"]
cand_timing_60 <- cand_timing_60[names(cand_timing_60) != "null"]

build_obj3_base <- function(cand_list){

  tab <- tibble::tibble(
    model_id  = names(cand_list),
    Deviance  = sapply(cand_list, stats::deviance),
    K         = sapply(cand_list, function(m) length(stats::coef(m))),
    AICc      = sapply(cand_list, AICcmodavg::AICc),
    R2        = sapply(cand_list, function(m) summary(m)$r.squared),
    `slope (SE)` = purrr::imap_chr(cand_list, ~{
      broom::tidy(.x) |>
        dplyr::filter(term != "(Intercept)") |>
        dplyr::mutate(txt = paste0(
          sprintf("%.3f", estimate),
          " (", sprintf("%.3f", std.error), ")"
        )) |>
        dplyr::pull(txt) |>
        paste(collapse = "; ")
    })
  ) |>
    dplyr::arrange(AICc) |>
    dplyr::mutate(
      `ΔAICc` = AICc - min(AICc),
      w = exp(-0.5 * `ΔAICc`) / sum(exp(-0.5 * `ΔAICc`)),
      Deviance = round(Deviance, 3),
      AICc = round(AICc, 2),
      `ΔAICc` = round(`ΔAICc`, 2),
      w = round(w, 3),
      R2 = round(R2, 3)
    ) |>
    dplyr::select(
      model_id, Deviance, K, AICc, `ΔAICc`, w, R2, `slope (SE)`
    )

  tab
}

model_key <- tibble::tribble(
  ~model_id,      ~Group,               ~Model,
  "meanT_30",     "Temperature",        "Mean temperature (30 days)",
  "maxT_30",      "Temperature",        "Maximum temperature (30 days)",
  "minT_30",      "Temperature",        "Minimum temperature (30 days)",

  "rainDays_30",  "Rain",               "Rainy days (30 days)",
  "totalR_30",    "Rain",               "Total rainfall (30 days)",

  "meanW_30",     "Wind",               "Mean wind speed (30 days)",
  "maxW_30",      "Wind",               "Maximum wind speed (30 days)",

  "meanT_60",     "Temperature",        "Mean temperature (60 days)",
  "maxT_60",      "Temperature",        "Maximum temperature (60 days)",
  "minT_60",      "Temperature",        "Minimum temperature (60 days)",

  "rainDays_60",  "Rain",               "Rainy days (60 days)",
  "totalR_60",    "Rain",               "Total rainfall (60 days)",

  "meanW_60",     "Wind",               "Mean wind speed (60 days)",
  "maxW_60",      "Wind",               "Maximum wind speed (60 days)",

  "lagTotR",      "Lagged effects",     "Previous year's total rainfall",
  "lagRainD",     "Lagged effects",     "Previous year's rainy days",
  "lagMeanT",     "Lagged effects",     "Previous year's mean temperature",

  "temp_rain",    "Joint effects",      "Mean temperature + rainy days",
  "temp_wind",    "Joint effects",      "Mean temperature + mean wind speed",
  "rain_wind",    "Joint effects",      "Rainy days + mean wind speed",

  "temp_x_rain",  "Interactive effects","Mean temperature × rainy days",
  "minT_x_wind",  "Interactive effects","Minimum temperature × mean wind speed",
  "lagR_x_temp",  "Interactive effects","Previous year's rainfall × mean temperature"
)

make_grouped_table <- function(cand_list, caption_text){

  group_order <- c(
    "Rain", "Temperature", "Wind",
    "Lagged effects", "Joint effects", "Interactive effects", "Other"
  )

  tab <- build_obj3_base(cand_list) |>
    dplyr::left_join(model_key, by = "model_id") |>
    dplyr::mutate(
      Group = dplyr::coalesce(Group, "Other"),
      Model = dplyr::coalesce(Model, model_id),
      Group = factor(Group, levels = group_order)
    ) |>
    dplyr::arrange(Group, AICc) |>
    dplyr::select(Group, Model, Deviance, K, AICc, `ΔAICc`, w, R2, `slope (SE)`)

  group_sizes <- tab |>
    dplyr::count(Group, name = "n") |>
    dplyr::filter(n > 0) |>
    dplyr::mutate(
      start = cumsum(dplyr::lag(n, default = 0)) + 1,
      end = cumsum(n)
    )

  kb <- knitr::kable(
    tab |> dplyr::select(-Group),
    caption = caption_text,
    align = c("l","r","r","r","r","r","r","r","l"),
 Twbooktabs = TRUE
  ) |>
    kableExtra::kable_styling(full_width = TRUE)

  for(i in seq_len(nrow(group_sizes))){
    kb <- kb |>
      kableExtra::pack_rows(
        as.character(group_sizes$Group[i]),
        group_sizes$start[i],
        group_sizes$end[i],
        bold = TRUE
      )
  }

  kb
}

make_grouped_table(
  cand_timing_30,
  "Table 3. Competitive candidate models (ΔAICc ≤ 4) explaining laying date (day of year) at Dronfield using 30-day pre-laying weather cues. K = number of parameters; AICc = small-sample AIC; ΔAICc = difference from the best model; w = Akaike weight. R² is shown for all models; slope (SE) is reported for each weather predictor."
)
Table 3. Competitive candidate models (ΔAICc ≤ 4) explaining laying date (day of year) at Dronfield using 30-day pre-laying weather cues. K = number of parameters; AICc = small-sample AIC; ΔAICc = difference from the best model; w = Akaike weight. R² is shown for all models; slope (SE) is reported for each weather predictor.
Model Deviance K AICc ΔAICc w R2 slope (SE)
Rain
Rainy days (30 days) 404580.3 2 12360.96 1024.96 0.000 0.161 -2.103 (0.126)
Total rainfall (30 days) 415698.5 2 12400.54 1064.54 0.000 0.138 -0.286 (0.019)
Temperature
Minimum temperature (30 days) 230069.4 2 11529.94 193.94 0.000 0.499 -5.010 (0.131)
Mean temperature (30 days) 241441.3 2 11600.33 264.33 0.000 0.475 -5.478 (0.151)
Maximum temperature (30 days) 280751.2 2 11820.41 484.41 0.000 0.389 -4.775 (0.157)
Wind
Maximum wind speed (30 days) 449433.4 2 12506.89 1170.89 0.000 0.022 4.154 (0.727)
Mean wind speed (30 days) 452345.2 2 12516.32 1180.31 0.000 0.016 4.397 (0.915)
Lagged effects
Previous year's total rainfall 444303.8 2 12278.29 942.29 0.000 0.010 0.013 (0.003)
Previous year's rainy days 447834.4 2 12289.62 953.61 0.000 0.003 0.062 (0.032)
Previous year's mean temperature 448801.1 2 12292.70 956.70 0.000 0.000 -0.456 (0.576)
Joint effects
Mean temperature + rainy days 203046.6 3 11349.65 13.65 0.001 0.558 -5.016 (0.141); -1.508 (0.091)
Mean temperature + mean wind speed 236458.4 3 11571.91 235.91 0.000 0.485 -5.453 (0.150); 3.669 (0.662)
Rainy days + mean wind speed 369830.8 3 12224.48 888.48 0.000 0.195 -2.169 (0.120); 4.930 (0.828)
Interactive effects
Mean temperature × rainy days 200878.1 4 11336.00 0.00 0.999 0.563 -4.428 (0.204); 0.693 (0.563); -0.159 (0.040)
Previous year's rainfall × mean temperature 235450.7 4 11373.62 37.62 0.000 0.476 0.026 (0.015); -4.875 (0.495); -0.001 (0.001)
Minimum temperature × mean wind speed 215499.8 4 11438.51 102.51 0.000 0.531 -1.416 (0.716); 11.636 (1.371); -1.198 (0.232)
make_grouped_table(
  cand_timing_60,
  "Table S3. Competitive candidate models (ΔAICc ≤ 4) explaining laying date (day of year) at Dronfield using 60-day pre-laying weather cues. K = number of parameters; AICc = small-sample AIC; ΔAICc = difference from the best model; w = Akaike weight. R² is shown for all models; slope (SE) is reported for each weather predictor."
)
Table S3. Competitive candidate models (ΔAICc ≤ 4) explaining laying date (day of year) at Dronfield using 60-day pre-laying weather cues. K = number of parameters; AICc = small-sample AIC; ΔAICc = difference from the best model; w = Akaike weight. R² is shown for all models; slope (SE) is reported for each weather predictor.
Model Deviance K AICc ΔAICc w R2 slope (SE)
Rain
Rainy days (60 days) 365166.2 2 12211.31 1450.61 0.000 0.243 -1.350 (0.062)
Total rainfall (60 days) 383366.2 2 12282.33 1521.62 0.000 0.205 -0.160 (0.008)
Temperature
Minimum temperature (60 days) 158946.3 2 10990.38 229.67 0.000 0.654 -5.499 (0.105)
Mean temperature (60 days) 172835.2 2 11112.60 351.89 0.000 0.624 -6.105 (0.124)
Maximum temperature (60 days) 209683.4 2 11394.57 633.86 0.000 0.544 -5.834 (0.140)
Wind
Maximum wind speed (60 days) 457772.5 2 12533.72 1773.01 0.000 0.004 -1.889 (0.803)
Mean wind speed (60 days) 459407.0 2 12538.92 1778.21 0.000 0.000 -0.572 (0.998)
Lagged effects
Previous year's total rainfall 444303.8 2 12278.29 1517.58 0.000 0.010 0.013 (0.003)
Previous year's rainy days 447834.4 2 12289.62 1528.91 0.000 0.003 0.062 (0.032)
Previous year's mean temperature 448801.1 2 12292.70 1531.99 0.000 0.000 -0.456 (0.576)
Joint effects
Mean temperature + rainy days 136023.4 3 10765.17 4.46 0.097 0.704 -5.383 (0.116); -0.798 (0.040)
Mean temperature + mean wind speed 171134.4 3 11100.18 339.48 0.000 0.628 -6.150 (0.124); 2.328 (0.612)
Rainy days + mean wind speed 336138.6 3 12085.12 1324.41 0.000 0.268 -1.389 (0.060); -1.748 (0.855)
Interactive effects
Mean temperature × rainy days 135421.3 4 10760.71 0.00 0.903 0.705 -5.851 (0.217); -1.503 (0.280); 0.045 (0.018)
Previous year's rainfall × mean temperature 167813.3 4 10889.01 128.31 0.000 0.626 0.067 (0.013); -4.329 (0.379); -0.004 (0.001)
Minimum temperature × mean wind speed 153677.5 4 10945.22 184.51 0.000 0.666 -2.223 (0.673); 11.131 (1.744); -1.092 (0.218)
ggplot(timing_weather_60, aes(x = lay_DOY)) +
  geom_histogram(bins = 15, fill = "grey70", colour = "black") +
  labs(
    x = "Laying date (day of year)",
    y = "Frequency",
    title = "Distribution of laying dates"
  ) +
  theme_bw() +
  theme(panel.grid = element_blank())

timing_weather_30_filt <- timing_weather_30 |>
  dplyr::filter(year >= 1993, year <= 2024)

ggplot(timing_weather_30_filt, aes(x = year, y = lay_DOY)) +
  geom_point(
    size = 1.8,
    colour = "grey30",
    alpha = 0.6
  ) +
  geom_smooth(
    method = "lm",
    se = TRUE,
    colour = "grey",
    linewidth = 1
  ) +
  scale_x_continuous(
    breaks = seq(1995, 2025, by = 5),
    limits = c(1993, 2024)
  ) +
  labs(
    x = "Year",
    y = "Laying date (day of year)",
    title = "Annual variation in laying date at Dronfield"
  ) +
  theme_bw() +
  theme(
    panel.grid = element_blank()
  )+
  geom_jitter(
  width = 0.2,
  height = 0,
  size = 1.6,
  colour = "grey30",
  alpha = 0.6
)
`geom_smooth()` using formula = 'y ~ x'
Warning: Removed 56 rows containing missing values or values outside the scale range
(`geom_point()`).

# Figure 3.3: Temperature × rainfall interaction with prediction ribbons

m_int <- cand_timing_30$temp_x_rain

df <- timing_weather_30 |>
  dplyr::filter(year >= 1993, year <= 2024)

# rainfall quantiles
rain_q <- quantile(df$rain_days_30,
                   probs = c(0.25, 0.5, 0.75),
                   na.rm = TRUE)

# prediction grid
newdat <- expand.grid(
  mean_temp_30 = seq(
    min(df$mean_temp_30, na.rm = TRUE),
    max(df$mean_temp_30, na.rm = TRUE),
    length.out = 80
  ),
  rain_days_30 = as.numeric(rain_q)
)

# predictions + SE
pred <- predict(m_int, newdata = newdat, se.fit = TRUE)

newdat$fit <- pred$fit
newdat$lo  <- pred$fit - 1.96 * pred$se.fit
newdat$hi  <- pred$fit + 1.96 * pred$se.fit

newdat$Rain_group <- factor(
  newdat$rain_days_30,
  labels = c("Low rainfall (25%)",
             "Medium rainfall (50%)",
             "High rainfall (75%)")
)

# plot
ggplot(df, aes(x = mean_temp_30, y = lay_DOY)) +
  geom_point(
    colour = "grey75",
    size = 1.6,
    alpha = 0.35
  ) +
  geom_ribbon(
    data = newdat,
    aes(
      x = mean_temp_30,
      ymin = lo,
      ymax = hi,
      group = Rain_group
    ),
    alpha = 0.12,
    inherit.aes = FALSE
  ) +
  geom_line(
    data = newdat,
    aes(
      x = mean_temp_30,
      y = fit,
      linetype = Rain_group
    ),
    linewidth = 1.2,
    inherit.aes = FALSE
  ) +
  labs(
    x = "Mean temperature (30 days pre-laying)",
    y = "Laying date (day of year)",
    linetype = "Rainfall conditions",
    title = "Temperature–rainfall interaction predicts laying date"
  ) +
  theme_bw() +
  theme(
    panel.grid = element_blank(),
    legend.position = "right"
  )

Laying date advanced with increasing pre-laying temperature, but this effect was contingent on rainfall conditions, indicating that thermal cues did not operate independently of broader environmental context. Warmer conditions are likely to reduce thermoregulatory costs and facilitate physiological readiness for breeding; however, the strength of this temperature effect varied with rainfall. Rather than rainfall acting as a direct trigger for breeding, wetter pre-laying conditions may reflect short-term improvements in ecosystem productivity or adult condition, thereby modifying how effectively temperature cues translate into earlier breeding.